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Abstract 

This  thesis  examines  a  remote  sensing  technique  for  measuring  the  atmospheric  struc¬ 
ture  constant  (C,^)  as  a  function  of  altitude  by  performing  spatial  correlation  of  wave  front 
sensor  measurement  5.  Two  point  sources  are  used  to  irradiate  two  wave  front  sensors  in 
the  aperture  plane  of  an  optical  system.  The  geometric  relationship  between  the  sources 
and  the  sensors  gives  rise  to  crossed  optical  paths.  At  the  point  where  the  paths  cross, 
the  correlation  value  of  the  turbulence  contribi'*ions  will  be  at  a  peak.  The  correlation  is 
shown  to  be  mathematically  related  to  the  structure  constant  in  terms  of  an  integral  of 
multiplied  by  a  path  weighting  function.  It  is  shown  that  the  path  weighting  function  ran 
be  made  to  have  the  characteristics  of  a  sampling  function  and  the  value  of  the  structure 
constant  can  be  directly  inferred  from  the  correlation  measurement.  The  vertical  resolution 
and  signal-to-noise  ratio  are  calculated  for  a  sample  case  of  two-layer  turbulence. 


Sensing  Refractive  Turbulence  Profiles  Using  Wave  Front  Slope  Measurements 

From  Two  Reference  Sources 


I.  Introduction 


1.1  Atmospheric  Turbulence 

An  electromagnetic  wave  traveling  through  the  Earths’  atmosphere  is  subject  to  a 
phenomena  known  as  atmospheric  turbulence.  Atmospheric  turbulence,  caused  by  differ¬ 
ential  heating  of  the  Earths’  surface  by  the  sun,  is  characterized  by  random  variations 
in  atmospheric  temperature  and  pressure.  An  electromagnetic  wave  passing  through  tur¬ 
bulence  experiences  a  bending  or  distortion  of  the  wave  front  shape.  Optical  systems 
operating  in  the  presence  of  the  atmosphere  are  subject  to  the  effects  of  the  turbulence. 
The  distortion  reduces  the  resolution  the  astronomer  achieves  when  looking  through  a 
telescope  and  disrupts  optical  communications.  The  familiar  twinkling  of  starlight  is  an 
example  of  the  distortion  caused  by  the  atmosphere.  The  light  waves  are  being  bent  as 
they  pass  through  the  atmosphere  to  the  Earth. 

The  efforts  to  mitigate  the  effect  of  atmospheric  turbulence  on  optical  systems 
spawned  a  branch  of  optics  known  as  adaptive  optics.  In  adaptive  optics,  techniques 
are  employed  to  reverse  the  effects  of  the  atmosphere  on  the  optical  wave  (7).  One  com¬ 
mon  technique  is  to  use  a  deformable  mirror.  A  control  system  senses  the  deformation  of 
the  incoming  wave  front  and  sends  signals  to  the  mirror.  The  mirror  surface  deforms  to 
a  shape  that  is  the  inverse  of  the  distortion.  When  the  distorted  wave  hits  the  deformed 
mirror,  the  distortion  is  nullified. 

Inherent  in  the  techniques  of  adaptive  optics  is  the  need  to  be  able  to  characterize  the 
nature  of  the  wave  front  distortion.  To  characterize  the  wave  front  distortion,  it  is  necessary 
to  know  how  the  turbulence  is  distributed  along  the  optical  path.  A  parameter  known  as 
the  atmospheric  structure  constant  is  commonly  used  as  a  measure  of  the  strength  of  the 
turbulence.  The  structure  constant  is  introduced  in  the  following  section. 
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1.2  Characlcnzulion  of  Turbulence 

The  atmosphere  may  l)e  thoiigiit  of  as  a  random  field.  'J'iie  pressure  and  lemj)cra- 
tnre  change  as  a  function  of  position.  Therefore,  the  pre.ssurc  and  temperature  are  non 
uniformly  distributed  along  an  optical  propagation  path.  In  the  theory  of  turbidence,  the 
strength  of  the  turbulence  is  characterized  by  the  magnitude  of  the  power  spectral  den 
sity  (psd)  of  the  index  of  refraction,  n.  This  psd  is  designated  as  *l’„(A’)  where  A'  is  the 
wavenumber  vector. 

The  magnitude  of  ‘h,i(A')  is  proportional  to  a  parameter  known  as  the  atmospheric 
structure  constant  of  the  refractive  index  fluctuations.  The  atmospheric  structure  constant 
is  denoted  as  C*  and  is  a  function  of  position  along  the  optical  path.  .Since  the  magnitude 
of  the  psd  is  proportional  to  C*.  C*  is  used  ;is  a  measure  of  the  strength  of  the  turbulence. 

1.3  Measurement  of  C" 

There  have  been  numerous  methods  used  to  mca.>ure  C*  (3.  10. 12,  M,  17).  .Actually, 
the  use  of  the  word  moiisuremcnt  is  a  misnomer  since  the  structure  constant  can  not  be 
directly  ntcastired  with  any  kind  of  instrument.  The  technique  u.scd  to  evaluate  C,*  is  to 
measure  a  parameter  of  an  optical  system  operating  in  the  presence  of  turbulence  and 
extr.act  the  value  of  C,®  from  the  measured  quantity.  Tried  (3)  was  one  of  the  first  to  apply 
this  approach  b\  performing  aspatial  cro.ss  correlation  of  stellar  .scintillation  measurements. 
As  with  Fried's  work,  current  techniques  rely  primarily  on  the  spatial  or  temporal  cro.sb 
correlation  |>roj)erties  of  the  intensity  of  the  o|)ticai  field.  This  research  will  lake  a  new 
approach,  based  not  on  the  intensity  of  the  optical  field,  but  on  the  correlation  properties 
of  the  wave  front  pha.se. 

l.f  .Approach 

Co.'isidcr  two  point  .sources,  denoted  pi  and  pn.  which  are  .separated  by  a  di.staiire 
A/j  as  shown  in  figure  l.I.  Two  wa\e  front  .beicsors  arc  placed  in  the  aperture  plane  of  the 
optical  .systeii!.  The  geometric  relationship  between  the  .sources  and  the  wave  front  .seii.iors 
gives  rise  to  cto.s.scd  optical  paths.  The  light  from  each  point  source  travels  a  dilferent 
optical  path  to  the  wave  front  sensors.  .\t  the  point  where  the  path.s  cro.ss.  the  turbulence 
contributions  from  the  two  souices  will  be  highly  correlated.  'J'liis  correlation  i.s  exploited 
to  calcidalc  C'~.  A  three  step  |>roce.ss  is  used  in  this  lhesi.s.  It  is  first  shown  that  the 
structure  constant  is  related  to  the  corielation  of  the  measurements  from  the  two  wave 


ilCf 


front  sensors  by 

j  dzCl{z')v:{:')  (1.1) 

where  C,  is  the  correlation  of  wave  front  sensor  measurement  i  with  wave  front  sensor 
measurement  2.  s'  is  a  position  along  the  vertical  path,  and  w{sf)  is  a  path  wdghting 
.function.  The  next  step  proves  that  the  weighting  function  can  be  shown  to  approximate 
a  Di'ac  delta  function  under  certain  conditions.  The  sifting  property  of  the  Dirac  delta 
function  may  be  used  to  evaluate  C*  from  the  relationship  in  equation  1.1.  Finally,  the 
signal-tc-noisc  ratio  (SNR)  of  the  correlation  measurement  is  calculalc^d  to  evaluate  the 
accuracy  of  the  method. 

1.5  Overview 

Chapter  II  reviews  the  background  of  atmospheric  turbulence  sensing  as  it  applies  to 
this  rcsearclu  Chapter  III  explains  the  mclhodoiegy  employed  in  this  rcscarcli.  Chapter 
IV  is  a  presentation  and  analysis  of  the  results.  Chapter  V  contains  concliisioris  and 
recommendations  for  further  study. 
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II.  Background 


2.1  Introduction 

This  chapter  examines  the  theory  of  atmospheric  turbulence  and  reviews  a  method 
used  by  Pried  to  calculate  the  strength  of  the  turbulence.  The  crossed-beam  method  of 
remote  sensing,  and  its’  application  to  this  research,  is  discussed. 

2.2  Theory 

2.2.1  Atmospheric  Turbulence  The  index  of  refraction,  n,  for  a  medium  is  the  ratio 
of  the  speed  of  light  in  a  vacuum  to  that  in  the  medium.  In  the  Earth’s  atmosphere,  n  is 
a  function  of  pressure,  temperature,  and  humidity.  As  such,  it  varies  with  height  above 
the  Earth’s  surface  and  n  may  be  expressed  as  a  function  of  position.  The  variation  of  the 
index  of  refraction  is  a  random  quantity.  The  atmosphere  may  be  visualized  as  consisting 
of  randomly  sized  pockets  of  air  called  eddies.  Typically,  eddies  range  in  size  from  a  few 
millimeters  to  tens  of  meters.  Each  eddie  is  at  a  different  pressure  and  temperature  and 
thus  has  a  different  n.  As  an  optical  wave  front  propagates  through  the  atmosphere,  the 
variation  in  n  from  eddie  to  eddie  distorts  the  wave  front. 

Since  the  characteristics  of  the  eddies  are  random,  n  is  a  temporal  and  spatial  random 
process.  As  a  result,  it  is  necessary  to  use  statistical  methods  to  describe  the  characteristics 
of  n.  In  atmospheric  optics,  it  is  common  to  characterize  the  second  order  statistics  of 
n  using  the  structure  function.  Using  angle  brackets  to  denote  the  ensemble  average,  the 
structure  function  is  given  by  (8:526-527) 

D„if^  =<|  n(f -b  ri)  -  n(fi)  |^>  (2.1) 


where  r  is  a  position  vector. 

2.2.2  Power  Spectral  Density  The  power  spectral  density  (psd)  of  the  index  of 
refraction  fluctuations,  is  a  measure  of  the  strength  of  the  turbulence.  As  Goodman 

states,  (6:388)  the  psd  may  be  regarded  as  a  measure  of  the  relative  occurrence  of  eddies 
with  dimensions  T  —  7?  is  called  the  wavenumber  vector  and  may  be  interpreted  as 

the  number  of  eddies  that  occur  in  a  unit  length.  As  Tatarski  shows  (15:19-21),  the  psd 
may  be  related  to  the  structure  function  by 

D„{iO  =2  J  J  J"jl  -  cos(7?  •  r)]  $„(/?)  d^K  (2.2) 
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2.2.3  Kolmogoroff  Spectrum  In  a  classic  paper  on  statistical  theory  (9),  Kolmogo- 
roff  developed  a  theory  of  turbulence  that  has  been  used  to  model  the  effects  of  the  atmo¬ 
sphere.  In  this  theory,  the  eddies  are  characterized  by  two  size  parameters;  the  outer  scale 
Lo  and  the  inner  scale  /<,.  If  the  eddie  size  is  between  Lo  and  /<,  (i.e.  ^  <  K  <  ^)i  the 
turbulence  is  essentially  isotropic  and  this  region  is  known  as  the  inertial  subrange.  Since 
the  turbulence  is  isotropic,  is  only  a  function  of  the  magnitude  of  K  and  the  spectrum 
is  expressed  as  (6:386-390) 

$„(/0  =  0.033  (2.3) 

where  K  =|  K  \.  The  quantity  is  known  as  the  structure  constant  of  the  refractive 
index  fluctuations  and  is  used  as  a  measure  of  the  strer  .  h  of  the  turbulence. 

2.3  Evaluation  of 

It  is  not  possible  to  measure  directly,  however,  the  evaluation  of  as  a  function  of 
position  has  been  widely  studied.  Both  vertical  and  horizontal  profiles  have  been  obtained, 
however,  the  vertical  profile  is  of  primary  interest  for  celestial  imaging.  Figure  2.1,  adapted 
from  tabular  data  in  reference  (1:2179),  is  a  plot  of  a  typical  vertical  profile  of  C*  from  1 
to  20  kilometers. 


Figure  2.1.  C,'  As  A  Function  Of  Altitude 
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In  1968,  PeskofT  (13)  published  a  mathematical  basis  for  calculating  C-l-  He  pro¬ 
posed  gathering  irradiance  data  from  optical  scintillations  and  correlating  the  irradiance 
me<rsurements.  The  correlation  value  is  used  to  extract  C,^.  At  the  same  time,  Fried  pro¬ 
posed  a  similar  experimental  implementation  (3).  Fried's  work  is  a  good  example  of  the 
correlation  method  and  is  explained  in  the  following  section. 


2.3.1  Correlation  of  Irradiance  Mcasui'enients  Fried  proposed  using  two  telescopes, 
each  having  a  restricted  field  of  view,  to  track  scintillating  stars  and  gather  irradiance  data 
(3:416-418).  The  intensities  measured  by  the  telescopes  are  denoted  /i(/.)  and  i->{t)  where 
i  is  time.  The  logarithmic  amplitude,  /,(/.),  is  computed  from  the  relationship 


b(0 

.<  u(0  >. 


(2.4) 


where  the  angle  brackets  denote  a  time  average.  In  turn,  the  spatial-temporal  log  amplitude 
covariance  C((/?,  r)  is  calculated  from  the  log-amplitudes  by 


-<  (/i(0“  ^i(0  >]  [^2(f  +  ^)~  <  ^2(0  >]  >  (2-5) 


Tatarski’s  work  (15:110-113)  contains  an  equation  which  relates  Ct{p.T)  to  C'^.  Fried 
applies  this  equation  and  reduces  the  resultant  integral  to  the  following  form: 


(2.6) 


In  this  expression,  F’(^)  is  a  function  that  represents  a  grouping  of  all  the  mathematical 
terms  that  depend  on  This  allows  C'^  and  2:  to  stand  alone  in  the  integral.  Fried 
analytically  .solves  the  function  /'( and  provides  a  table  of  values.  Using  this  table,  the 
measured  Ci{p)  is  now  related  to  the  unknown  by  a  linear  integral  equation. 

The  approach  used  by  Fried  has  become  .somewhat  standard.  .Measuiemcnts  are  made 
of  a  property  of  an  optical  field  and  the  measured  {|uantities  are  lelatcd  to  C',-|  (usually  In 
an  integral).  The  terms  in  the  equation  are  grouped  into  a  separate  function  to  allow 
to  stand  alone.  'I'he  &ei)arate  function  is  called  a  path  weighting  function  and  relates  C',‘ 
to  the  measured  quantity. 


2.3.2  The  Crossed- lk:am  Method  In  1967,  Fisher  and  Krau.se  (2)  propo.scd  a  method 
tocxtract  turbulence  values  from  measurements  made  using  rros.sed  laser  beams.  'J'no  laser 


beams  intersect  at  a  right  angle  in  a  turbulent  medium.  At  the  point  of  intersection,  there 
is  a  related  (correlated)  fluctuation  in  both  beams.  Away  from  the  intersection  point, 
the  turbulence  effects  are  uncorrelated  and  produce  uncorrelated  effects  on  the  beams.  If 
the  covariance  is  calculated,  the  uncorrelated  effects  yield  a  zero  average  value  while  the 
correlated  effects  yield  a  non-zero  average  value.  The  correlation  value  at  the  intersection 
is  a  function  of  the  magnitude  of  the  turbulence  fluctuations  at  or  near  the  intersection 
point.  The  covariance  values  can  be  used  to  evaluate  the  strength  of  the  turbulence.  They 
show  the  method  produces  accurate  results  for  measuring  turbulence  in  the  shear  layer  of 
a  subsonic  jet. 

In  1974,  Wang  and  his  colleagues  (17)  applied  the  crossed  beam  theory  to  sense 
horizontal  wind  and  refractive  turbulence  profiles.  They  did  not  require  the  beams  to  be 
perpendicular.  Calculating  the  covariance  function  of  intensity  fluctuations  of  the  crossed- 
path  signals,  their  development  yields  an  expression  very  similar  to  equation  1.1.  is 
related  to  the  covariance  by  an  integral  involving  a  weighting  function.  They  profile  the 
weighting  functions  for  three  different  transmitter-receiver  configurations. 

2.^  Conclusion 

Characterizing  the  strength  of  atmospheric  turbulence  by  remote  sensing  has  been 
widely  studied.  The  crossed-beam  method  is  often  used  to  measure  a  parameter  of  an 
optical  field.  The  measured  values  are  related  to  the  atmospheric  structure  constant,  and 
the  structure  constant  serves  as  a  me^lsure  of  the  strength  of  the  turbulence.  Much  of  the 
]L..ot  work  has  concentrated  on  optical  scintillation  measurements  wrere  the  intensity  of 
the  light  is  the  quantity  of  interest.  It  is  proposed  that  ciossed  beams  may  be  used  to 
obtain  a  measurement  of  wave  front  phase,  rather  than  intensity,  and  the  phase  correlation 
properties  may  be  used  to  calculate  C,^.  The  next  chapter  presents  the  development  of  this 
proposal. 
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III.  Methodology 


3.1  Introduction 

As  discussed  in  the  previous  chapter,  a  primary  method  used  to  obtain  vertical 
profiles  of  C,^  has  been  through  the  use  of  optica!  scintillation  measurements.  These  mea¬ 
surements  utilize  the  intensity  of  the  optical  field.  Instead  of  using  intensity  as  the  quantity 
of  interest,  it  is  shown  in  this  chapter  that  measurements  of  wave  front  phase  slope  may 
be  used  to  obtain  profiles  of  C^. 

Consider  two  point  sources,  located  at  points  pi  and  p2)  <ind  separated  by  a  vector 
distance  A/j,  as  shown  in  figure  3.1.  Two  wave  front  sensors,  and  are  placed 


4- 


in  the  aperture  plane  of  the  optical  system.  The  wave  front  sensors  are  separated  b\  the 
vector  distance  Aa:.  Each  sensor  measures  the  slope  of  the  distorted  wave  front  across  a 
finite  area  called  a  subaperture.  Each  sensor  is  diiided  into  multiple siibaperturcs,  each  of 
which  provides  a  slope  measurement. 


3-1 


It  is  common  to  normalize  the  aperture  weighting  functions  (16:1771)  so  that 

jd-xWn{x)  =  l  (3.1) 

where  n,  in  the  case  of  figure  3.1,  may  be  1  or  2. 

The  point  sources  are  at  a  height  z,  above  the  aperture  plane  and  z'  is  a  point  along 
the  z  axis.  As  can  be  seen  from  the  figure,  the  geometric  relationship  between  the  sources 
and  the  sensors  gives  rise  to  crossed  optical  paths.  The  light  from  each  point  source  travels 
a  different  optical  path  to  the  wave  front  sensor.  At  the  point  where  the  paths  cross,  the 
turbulence  contributions  will  be  highly  correlated.  This  correlation  is  exploited  to  calculate 

Cl 

This  chapter  contains  mathematical  derivations  based  on  the  correlation  of  wave  front 
phase  slopes  for  the  geometry  of  figure  3.1.  The  first  derivation  proves  that  the  correlation 
of  wave  front  slopes,  denoted  C, ,  can  be  related  to  by  the  following  equation. 

C,  -  J  dz'  Clz')  w{z')  (3.2) 

where  z'  is  a  position  along  the  optical  path  and  w{z')  is  a  path  weighting  function.  If  w(z') 
has  the  shape  of  a  Dirac  delta  function  (i.e.  a  sampling  function),  the  sifting  property  may 
be  applied  to  equation  3.2  to  directly  relate  C,  to  Cl 

Following  the  derivation  of  equation  3.2,  the  expression  for  w{z')  is  simplified  so  that 
it  depends  only  on  terms  of  the  measurement  geometry  {Ax,  Ap,  z' ,  etc.).  This  allows 
v}{z')  to  be  plotted  to  examine  the  sampling  characteristics.  It  is  shov/n  in  Chapter  4  that 
the  path  weighting  function  derived  in  this  chapter  does  not  possess  an  ideal  sampling 
shape.  Therefore,  a  discussion  is  included  in  this  chapter  about  a  method  to  improve  the 
shape  of  w{z'). 

Finally,  an  expression  for  the  signal-to-noise  ratio  of  the  correlation  measurement  is 
derived.  This  is  used  to  quantify  the  accuracy  of  the  measurement  technique. 
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3.2  Overview  oj  Derivation 


Since  the  derivation  of  equation  3.2  is  rather  involved,  it  is  useful  to  first  outline  the 
procedure.  The  first  step  is  to  define  what  is  meant  by  the  correlation  of  wave  front  slopes. 
The  ne.xt  step  is  to  define  an  expression  which  relates  the  slope  signals  to  the  measurement 
geometry  of  figure  3.1.  This  allows  the  correlation  of  the  wave  front  slopes  to  be  directly 
related  to  the  measurement  geometry. 

After  the  expression  for  the  correlation  of  the  wave  front  slopes  is  developed,  the 
ensemble  average  is  taken.  This  is  done  because  the  slope  signals  are  defined  in  terms  of 
the  wave  front  phase  of  the  optical  wave  Irom  each  source,  which  is  a  random  process. 
Taking  the  ensemble  average  introduces  the  second  order  moment  of  the  wave  front  phmse. 
An  expression  in  the  literature  relates  the  second  order  moment  of  the  phase  to  the  phase 
structure  function.  Recalling  from  Chapter  2  that  the  structure  function  is  related  lo  C,^, 
the  next  step  is  to  use  an  e.xpression  (also  from  the  literature)  to  relate  the  two  quantities. 
The  resulting  expression  will  be  an  integral  over  z'  that  includes  and  several  other  terms. 
The  t^erms  other  than  Cl  are  grouped  and  called  «;(*'),  the  path  weighting  function. 

The  final  part  of  the  derivation  concentrates  on  the  path  weighting  function.  By 
assuming  Kolmogoroff  statistics  for  the  atmosphere,  is  reduced  lo  a  form  which  is 
only  dependent  upon  known  (pianlities  from  the  measurement  geometry.  It  is  therefore 
p0!3.sib!e  to  plot  w(z')  lo  e.xamine  how  closely  it  appio.ximates  a  sampling  function. 


3.8  Derivation  of  the  Path  Wcighlintj  Fvnetion 

This  section  contains  the  dcrivaiion  of  the  path  wnighting  function  denoted  as 
in  equation  3.2. 


3.. 3.1 


Correlation  of  .Slope  .Sifjnal.’t 


Tim  {I..-'t  step  is 


'.0  d'T:iie  the  correlation  as 


$■,) 


(3.3) 


where  .S]  and  .S2  arc  the  .slope  signals  from  two  Indieidual  wave  fionl  .scii.sor  sub<ipcrtuio.s 
separated  by  A.r:.  Recalling  the  goometr.',  of  figure  3.1 ,  .Sj  is  the  .slope  measurement  arising 
from  reference  source  I,  and  .s-.-  is  the  slope  measun'inenl  arising  from  reference  source  2. 
For  the  remainder  of  this  development,  it  is  .assumed  I  hat  each  wave  front  sensor  consists 
of  only  a  single  subaperture  and  the  .separation  between  the  sensors  is  A;;;. 
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The  next  step  is  to  find  an  expression  for  the  sloi  :  signal.  Referring  to  the  literature, 
(16:1771)  the  phase  of  the  incoming  wave  is  designated  where  z  is  a  two-dimensional 
position  vector  and  n  refers  to  the  source  of  the  wave  front.  For  the  geometry  of  figure  3.1, 
n  may  be  1  or  2.  It  is  convenient  to  define  a  z-.  mean  phase  (4:1914)  related  to  i’niS)  by 

4>„{x)  =  ”  /  o  '-  IK,'  9)  -tpnix')  (3.4) 

The  output  of  a  wave  front  sensor  is  a  noisy  •  leasurement  of  the  average  slope  of  ^n(^0 
over  the  aperture  and  is  designated  ,s„.  The  c'l  s'-oii  for  s„  is 

s„=  j  c  X  Vl’’„(f)  (V.^„(®)  •  d„]  o.„  (3.5) 

where  Wn{x)  is  the  aperture  weightiut,  function,  '7^r.(ic)  is  the  spatial  gradient  of  the 
phase  4>n,  d„  is  a  unit  vector  in  the  direction  of  sensitivity  of  the  sensor,  and  a„  is  the 
slope  measurement  error.  The  measurement  error  is  attributed  to  photon  noise  in  the 
detection  process.  Integration  by  parts  yields 

Sn  =  -J  [VW„{x)  ■  4]  (l>„{x)  -1-  (3.6) 

To  simplify  v.otatioh,  let  VW^ff)  •  d„  =?  hF*(«)  where  gradient  of  W„(x)  in 

the  direction  of  the  n*'*  sensor.  Rewrite: 

s„  =  -  /  drx  l'F’(x)  (f>„{x)  -f  a„  (3.7) 

Assuming  the  subaperture  weighting  functions  are  identical  except  for  being  shifted  relative 
to  each  other,  it  is  possible  to  write  H'j(£)  =  l'F(®)  and  V72(f)  =  W{?i -  Ax)  where  l'F(.r) 
is  the  subaperture  weighting  function.  Using  these  representations  and  plugging  equation 
3.7  into  equation  3.3, 


C,(A®)  = 

d-x  W*{xj  ^i{x)  -p  ttj  I 
X  drx'  W"{x'  —  Ax)  4>2{x')  -1- 


(3.8) 
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j^.fter  muitiplying  aiul  regrouping,  C,{Ax)  becomes 


C,(a1c)  = 

j d-x  J  d-x'  V7'(.t)  W\x'  -  Ax)  [(^,(z) 
+  j d:-xW^{x)<j,i{x)  a, 

+  j  d’x'  V7'(x'  -  Ax)  4--<{x')  ft] 

4"  Cl’i  0'2 


(3.9) 


3.3.2  Ensemble  Average  of  Cf  As  mentioned  in  the  Overview  equation  3.9  contains 
random  quantities  and  it  is  tlierei’ore  necessary  to  use  statistical  methods  to  i-ioceed  with 
the  derivation,  (j)  is  a.  two  dimensional  random  process,  a  is  assumed  to  be  a.  zero  mean 
random  variable  independent  of  (f>.  Co.nsider  the  ensemble  average  of  C.,(<^u-),  denoted 
with  the  angle  brackets: 

{Cs{Ax))  =  (3.10) 

J  d-x  J  d-x'  1'7'X.i:)  W^x'  -  Ax)  {(f>i{x)  <f>2{x')) 

+  0 
+  0 

+(0102) 


.Simplifying, 


(C’,(Aa;))=  (3.1  n 

J  d'x  J  d-x'  U''''(;c)  H''(a'  -  Ax){(l>]{x)  (.h.,(x')) 

+(o  102) 


Consider  the  correlation  of  the  noi.se  functions,  (01O2),  in  equation  3.11.  Recall  that  o> 
and  O2  i'le  noi.se  from  the  wave  front  sensor  nica.surements.  In  eciuation  29  of  rcfcrei.'e 
(16).  Wallncr  expresses  the  correlation  of  slope  measurement  noises  as 


where  n  and  n'  designate  the  wave  front  sensor  subapertures,  ct“  is  the  photon  noise 
density,  and  8{n-n')  is  the  Dirac  delta  function.  Equation  3.12  indicates  that  for  different 
subaperturcs  (i.e.  n  ^  n'),  (a„an')  =  0.  For  the  geometry  under  consideration,  the 
subapertures  do  not  coincide  and  the  noise  term  is  zero.  Wallner  also  derives  (16:1773- 
1774)  an  expression  for  the  second  order  moment  of  the  phase,  terms  of 

the  phase  structure  function  Dio: 

*')  +  (){x)  +  g{x')  -  a  (3.13) 


where 


and 


and 


Di2(.r,.'c0  =  {[i>{x)  -  tp{x')]-) 
g(x)  =  ^  J  "  i'F(a;")  D(x,x") 
a  =  ^  y  dx"  J  dx"'W{x")W{x"')D{x"  ,x^") 


(3.14) 

(3.15) 

(3.16) 


Since  W‘(x)  is  an  odd-symmetric  function,  the  integratin'  p  over  x"  and  x"'  yield  a  zero 
Vcdue.  Therefore,  all  terms  in  equation  3.13,  except  —  |Dii(x,x'),  will  go  to  zero. 

Using  equations  3.13  and  3.12  in  equation  3.11,  the  correlation  of  the  slope  signals 
may  be  written  as 


(C,(A'x))  =  (3.17) 

-^D,2(x,x') 

To  further  simplify  equation  3.17,  it  is  necessary  to  have  an  expression  for  Di2(x,x')  that 
is  applicable  to  atmospheric  turbulence  and  the  geometry  of  figure  3.1.  This  is  discussed 
in  the  next  section. 


J  drx  J  d‘x' 


ir  (x)  W\a7>  -  Ax) 
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3.3.3  Path  Weighting  Function  For  spatially  stationary,  transverse  isotropic  tur¬ 
bulence  for  the  two  source  geometry  of  figure  3.1,  Lutomirski  and  Buser  (11:2160)  show 
that  the  phase  structure  function  may  be  expressed  as 


Di-^x-x')  = 


Stt-V  r  dz'  rdJij.  H\  i\,ih\,0,z') 
Jo  Jo 

X  {l  -  Jo  (/a  |(Ap)  (^)  +  {x-  x')il  - 


(3.18) 


where  ^n{Kx,F\,z')  is  the  spectral  density  of  the  index  of  refraction,  A/i  is  the  vector 
separation  of  the  optical  sources,  and  x  and  x'  specify  two  points  in  the  aperture  plane  of 
the  optical  system.  The  spectral  density,  <1»„,  is  a  function  of  distance  from  the  aperture 
plane,  z'.  and  of  wavenumber  components  perpendicular  (7fj.)  and  parallel  (A',)  to  the 
axis.  Jo  is  the  zero-order  Bessel  function. 

Following  Welsh  (18).  assume  the  power  spectrum  of  the  index  of  refraction  is  sep¬ 
arable  in  the  perpendicular  and  parallel  wavenumber  components  (Kx,K.)  and  z'.  This 
yields 

=  <!>:,( A'a.,  A-,)  CUz')  (3.19) 

Substituting  equation  3.18  and  3.19  into  equation  3.17  yields 


(C,(A:r))  r. 

J  d-x  Id-?  lF-’(.i:)  H- ■'(:)?  ~  aIc) 

X  {-‘W-lr)  ["  dz'  r  dI{xKxKU<x.m1,(^') 
Jo  Jo 

)} 

'J'o  match  the  form  of  equation  3.2.  this  may  be  rewriltcn  as 


X  ^  I  —  ./n  J. 


(A/;)(-)-l-(.r-.f'){l--) 


(C,(A:r))=  r  dz'  C-[z')  w{z'.z,.^:v.,di,) 
Jo 


where  w{z',z,,^x,i\j))  is  a  path  weighting  function  and  is  defined  as 


(3.20) 


(3.21) 


(3.22) 


w{ 


■',2,,Ax,Ap)  = 

J  d-x  J  d-x'  l'P(x)  14^-' (aT'  -  ^x) 


X  {-4n-lc^)  dh\  R\  i>:.(A'x,0) 

X  (l  -  Jo  (a'x  |(A7a)  (^)  +  {x  -  x')il  -  ^)|)  } 


3.3.4  Evaluation  of  Path  Weighting  Function  To  further  evaluate  the  patli  weight¬ 
ing  function,  it  is  necessary  to  assume  a  form  for  $1,.  For  atmospheric  turbulence,  Kol- 
mogoroff  (9:154)  has  shown  that 

$;,(/ix,  J\\)  =  0.0337i  (3.23) 


where  K  =|  A'x  -f  Kt  |.  Using  3.23,  the  expression  for  the  path  weighting  function  may  be 
rewritten  as 


w{z\Zf,£x,Sj))  = 

J  (Ex  J  (Ex'  VF’(.'i;)  W-'(x'  -  Ax) 

X  r  (lli\  h\  0.033A"T^ 

Jo 

X  |l  -  Jo  (Ap)  (^)  -f-  (x  -  .•?)(]  -  7-)  )  I 


(3.24) 


Employing  the  integral  identity  (15:269) 


j(“  (fe(l  -  =  5  |■2'• 


.. .  “(7^-4) 
sin - - - 


-I 


where  1  <  p  <  3,  and  combining  terms, 


w{z',z,,Ax,Ap)  = 

-4(0.033)«  */;'  J  (Ex  J  (Ex'  W.’(:?  -  Ax) 


X  n 


2I  r--*(  )sin(^) 
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-1 


(^'70(r)  +  (*“  '’^'Ki  -  r) 


X 


Simplifying. 


w(2',  Zs ,  i^x,  Ap)  =  (3.2?) 

-i.m-  j  d-x  J  d-x'  W^{x)  W^{x'  -  Ax) 

3. 

X  (^7^  (^)  +  (x  -  .-cOCl  - 

To  evaluate  equation  3.27,  it  is  necessary  to  assume  an  e.xpression  for  IT’(f).  It  is  not 
uncommon  to  see  squar*.  .Mibapertures  in  a  wave  front  sensor  and  the  rectangle  function  is 
commonly  used  to  describe  the  tranoinittance  of  such  an  aperture  (5:67).  Mathematically, 
this  implies  that 

rect(y)  rect(y)  (3.2S) 

L“  Jj  L 

where  the  factor  is  nece.ssary  to  achieve  normalization  to  1  (see  equation  3.1).  Taking 
the  derivative  of  equation  3.2S  yields  H'’(a;),  as  required. 

In  general,  the  derivative  may  be  expre.ssed  using  the  delta  function  (5:63-66): 

VK*( .;)  =  VlV(.f) .  d  =  +-)-  S{x  -  |)1  rcct(|)  (3.29) 

Implicit  in  this  o.xprcssion  is  the  fact  that  the  slope  in  the  x  direction  is  the  quantity  under 
consideration  (i.e.  d  is  an  x  directed  unit  vector). 

S-9..1.J  Bvnlualion  of  Ahsotule  Vtdue  QmmlUy  It  is  convenient  at  this  point 
to  evaluate  the  absolute  value  quantity  in  equation  3.27.  Representing  the  vectors  by  their 
components  in  the  x  and  y  ilircctions 


Ap  =  Apri  +  Ap,jy, 

(3.30) 

Ax  =  -b  Ayy, 

(3.31) 

X  =  x.r  +  yy 

(3.32) 

x’  -  X  x-i-yy 

(3.33) 

and  substituting  into  equation  3.27  yields, 


{:{.:{■!} 


y-. y-. ^>)  = 

(Ap.i  +  +  \xx  +  yi)  -  {x'x  -f  y'y)l(l  - 

Simplifying. 

y(x.  X  j  T/.  t/  5  ,  i ,  =  (-j.-jo) 

+  (a:  -  a^'Kl  “  ^)1*  -r  t*  (^  “  ^  )( * 

A, 

Typically.  *,  is  iimch  greater  than  z'  over  the  range  of  interest.  'I’liis  is  easily  seen  if  the 
sources  are  laser  gnid^tars.  The  guidestar  altitude  would  he  in  the  range  of  90  to  100 
kilometers.  Recalling  figure  2.1.  the  range  of  interest  for  the  turbulence  is  well  below  this 
altitude.  Therefore,  r  is  much  less  than  1  and  it  is  reiusonablc  to  use  the  appro.xiination 
(1  —  1  and  simplify: 

/(l.  X  j  Jt;  y  .  ^Pxt  X  .  *»)  —  (3.-1G) 

4-  X  _  xf  +  (Ajvr  +  if  -  y')|' 

S.9.5  SimplificMion  of  Pnlh  Weighting  Function  It  is  now  possible  to  simplify  the 
path  weighting  function  in  equation  3.27.  Replacing  the  aperture  functions  with  thdr 
mathematical  representations. 

z,.  Sx.  £p)  =  (3.371 

-lAG!rJ.r-x  J,r-£'  f(x,^.y.i/,^p,.^p,.z^.z,)  x 

{  -b  ^  ^  llrcctC  1 1|  >: 

-)]  rcrtC— ■ 

where  /{x.x/.y.y'.Apr.^p^-z’.z,)  is  defined  by  ct]iiniioii  3.36.  Kniploying  the  .sifting 
property  of  the  Dirac  delta  function  allows  the  intograiioiis  «\er  j  .ind  s’  to  Ik*  perff»rm«I 
and  the  weighting  function  reduce  to 


3  in 


(3.38) 


w{z',z„Ax,Ap)  = 

-lA6k^L~‘*  J  dyx&cx(^)  J  dy'  x 

^,y,y\^Pr,^Py,z',^s)  - 

^,y,y'^^Px,^Py,z',z,)  + 
f{^,^x  +  ^,y,y'i^Px,^Py,z',z,) 

Let 

9{y  -  y')  =  (3.39) 

f{~,Ax-  -,y,y\Ap^,Apy,z\z,)- 

f{-^,^^-^^iy^y'APx,^Py,z\z,)- 

fi^>^^-^,y,y\^Px,i^Py,z',z>)-^ 

L  L 

fi~,Ax  +  -,y,  y\  Ajh,  Apy,z',  z,) 

where  dependence  of  g  on  L,  Apx,  Apy,  z',  and  z,  is  dropped  for  notational  convenience. 
Substituting  3.39  into  3.38  gives 

w(z',  z,,Ax,  Ap)  =  (3.40) 

-lAek~L~'^  j  dy  rectij)  j  dy' XGd,{—-^)  g{y  -  y') 

By  making  the  change  of  variables  u  =xy-y'  and  v  =  equation  3.40  can  be  rewritten 
as 

w{z',z,,Ax,Ap)  =  (3.41) 

-lAQk^L~‘^  J  du  J  dv  rect(-"2"^-)  rect(- — ^)Sf(u) 

Another  change  of  variables  results  in 
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(3.42) 


w(z',z,,Ax,Ap)  = 

-lAQk^L~'^  j  du  g(u)  j  dv'  rect(^)  rect(^^ - ^ — —) 

Evaluating  the  integral  over  v'  yields 

w{z',  z,,  Ax,  Ap)  —  -lA6k^L~^  J  du  g(u)  tri(^  (3.43) 

where  tri  is  the  triangle  function.  Using  equations  3.39  and  3.36  the  weighting  function 
can  be  written 

w(z',  Zg,  Ax,  Ap)  =  (3-44) 

-lA6k^L-^J  duirii^!^^)  X 

2  [(Ap^—  -  Ax)^  +  (Apy—  +  u)^ 

L  Zt  2$  i 

-  [(Apx—  -Ax-  Lf  +  (Apy—  +  uf 

-  -Ax  +  Lf  +  (Apy^  +  u)^] 

Without  loss  of  generality,  assume  the  apertures  and  point  source  geometry  is  such 
that  there  is  no  offset  in  the  y  direction.  This  implies  Ay  and  Apy  are  0  and  equation  3.44 
simplifies  to 

w{z',  z,.  Ax,  Ap)  =  (3.45) 

-lAGk^L~^  J dMtri(^)  x 

r  ,  £ 

2  {Aj)j- —  Ax)^  + 

,  ^9 

—  (Apc- —  Aa:  -  T)^  +  u^ 

2, 

z'  .  1  ^ 

-  (Apx - Aa:  +  L)'‘‘  + 

The  quantity  Ap^—  —  Ax  may  be  interpreted  as  the  vector  distance  between  the  crossed 
ray  paths  as  a  function  of  z'  (see  figure  3.1).  Let  p  =  Ap^^  -  Ax  and  rewrite  as 
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(3.46) 


Let  u'  =  This  allows  the  integral  to  become  dimensionless  and  the  integration  may  be 
performed  from  -1  to  1,  regardless  of  the  size  of  L.  To  simplify  notation,  the  weighting 
function  is  designated  as  ■w{z')  and  is  written  as: 


w{z')  = 


(3.47) 


-I  /"* 

L.46fc“L~  J  dll'  tri(w')  x 


2|(f)’  +  “"  -1(7-1)’  +  “'’  ■-1(^  +  1)’  +  “'’ 


This  is  the  final  e.xpression  for  the  path  weighting  function.  The  characteristics  of  this 
function  are  examined  in  the  next  chapter. 


SJ,  Improvement  of  the  Path  Weighting  Function 

It  is  shown  in  the  next  chapter  that  the  shape  of  the  path  weighting  function  derived 
in  the  previous  section  does  not  {)osscss  the  characteristics  of  a  sampling  function.  The  low 
spatial  frequencies  of  the  phase  perturbations  are  causing  a  common  tilt  in  both  slope  inca 
surements.  This  effect  shows  up  as  a  dc  offset  in  the  weighting  function.  The  method  used 
to  improve  the  shape  of  the  path  weighting  function  is  to  combine  the  measurements  for 
different  size  apertures.  This  action  should  remove  the  common  tilt  component.  Equation 
3.47  is  used  with  different  values  of  L  to  achieve  the  desired  result.  Plots  of  the  improved 
weighting  function  for  this  method  are  given  in  the  next  chapter. 


S.5  Signal-io- Noise  Ratio  of  the  Correlation  Measurement 

In  order  to  quantify  the  usefulness  of  the  techniques  developed  in  tliis  rc.search.  the 
signal-to- noise  ratio  (SNR)  of  the  correlation  measurement  C,  is  calculated.  In  general, 
the  SNR  may  be  expressed  as  (20:170-171) 
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where  the  quantity  in  the  denominator  is  recognized  as  the  square  root  of  the  variance  (i.e. 
the  standard  deviation).  The  derivation  of  the  variance  is  rather  lengthy,  however,  the 
procedure  is  very  similar  to  that  followed  in  the  beginning  of  this  chapter.  For  that  reason, 
the  derivation  of  the  expression  for  the  variance  is  summarized  in  the  following  section. 

3.5.1  Derivation  of  the  Variance  The  variance  of  C,  is  given  by  the  square  of  the 
denominator  of  equation  3.47 


VAR  =  {C^;iAx))  -  {C^iAx))^  (3.49) 

The  second  term  on  the  right  hand  side  of  equation  3.49  is  the  square  of  the  already 
derived  result  given  by  equations  3.21  and  3.47.  The  first  term  on  the  right  hand  side 
is  the  ensemble  average  of  the  square  of  C\.  To  calculate  this  first  term,  equation  3.9 
is  squared  and  the  ensemble  average  is  taken.  Squaring  equation  3.9  yields  16  terms. 
However,  a  is  a  zero  mean  random  process  and  12  terms  go  to  zero  when  the  ensemble 
average  is  taken.  Equation  3.49  may  thus  be  written  as: 


VAR  = 


I J 1 1  " 

X  kF‘'(f)  VF-X.-?)  -  Aa-)  -  Lx) 

X  {(p\{x)  (l)-.{x’)(f)a{x")<p.iix>")) 


(3.50) 


-// 

^// 


d'x  d'x'  M'’^(.i;)  W’{x')  (^i(.'e)  <l>2{x')){ai  aj) 

d~x"  d'x’"  lT''(;t"  -  A;i;)  -  A.t)  ( d)i{x")d>A{x"')){(.\-i  a-,) 


+  (<*!  ®i){f''2  ^*2) 


j  j  d-x£:?W''{x)  lF'(:i?)  {4>dx)<!>'.{P)) 


The  reader  will  notice  equation  3.50  is  similar  to  equation  3.10.  'J'he  fundamental  difference 
is  the  pre.sence  of  the  fourth  order  moment  of  the  phase  rather  than  the  .second  order 
moment.  Calculation  of  the  fourth  order  moment  of  the  pha,sc  is  possible  if  Gaussian 


statistics  are  assumed  since,  for  Gaussian  processes,  the  fourth  order  moment  may  be 
expressed  as  a  combination  of  the  first  and  second  order  moments  (6:39).  Using  the 
Gaussian  assumption,  the  variance  expression  takes  the  following  form: 


VAR  = 


2.13  X-®  r  r  dz'dz"  Cl{z')Cl{z") 

Jo  Jo 

X  J  du  tn(  j)  j  dv  tri(  j)  [/i(ii,  v,  z',  z")  +  2/2(11,  v,  z\  z")] 
+1.46  k'*  L~^  j  dz'  Cf^{z')  J  du  2/3(11,2') 


-  lAQ  L~^  J  dz' Cl{z')  j duin(^)fji{u^z') 


(3.51) 


The  functions  f\  to  /4  are  similar  in  form  to  g  of  equation  3.39.  Letting  p  equal  Apx^, 
the  functions  are  as  follows. 


fi{u,v,z',z")  = 


4  pp*-)  +  ^  J 

-2  [(Ap,^  +  Lf  +  11^]  ^  [(Ap.^)^  + 
-2  j^(Ap;c^)^  +  M‘j  ^Aj)x— -  L)- +  V- 
-2  (Ajj^^)- +  11-]  +  Lf  +  V- 

-2  {Apx— -  L)~  +  {Apx—)~  +  v- 


+  j^(A7J,^  +  L)' +  11'j  (Api^  +  L)-  +  v- 


.  r\2  .  .2 


+  [( A?;*  ~  +  L)-  +  11-J  |^( Aji^  ^  -  X)-  + 

+  [(Apr^  -  Lf  + 11- j  [(Ap*^  +  Lf  +  n-J 


(3.52) 
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f2iu,V,Z,Z  )  = 


4  |^(^Px—  - 


-2  [(Ap*-  -  A®  +.  Lf  +  ti-J  [(Apx^  -  Ax)2  + 

-2  j^(Api|-  -  Aa;)^  +  u^j  [^(Ap*^ —  Ax  +  Lf  + 

-2  [(Ap*—  -  Ax  -  Lf  +  i^Px- —  Ax)^  + 

l.  i  V 

-2  (Ap:c^  -  Ax)^  +  «*•  (Apx^ - Ax-i)^  +  x^ 

J  L 

+  [(Ap,:— -  Ax  +  i)^  +  w’  (Ap*^ Ax  +  i)2  +  t; 

+  [(Ap^— -  Ax  +  i)^  +  «^  (Apx— -  Ax  -  if +  x 

L  '  Xf  i  \.  z, 

+-[(Ap*^  -  Ax  -  Lf  +  ■  j^(Ap*^  -  Ax  +  i)®  +  V 

Ax  -  Lf  +  |^( 


fz{u,z')  = 

2f(Ap,-f +  „ 


(3.55) 


/,(«,  z'y= 


2 


(Apx - Arc  -  LY  + 

('Ap*—  -  Ao:  +  LY  + 


At  this  point,  in  order  to  get  some  simplified  results,  it  is  necessary  to  assume  a  form  for 
the  turbulence.  The  assumption  is  made  that  the  turbulence  is  confined  to  two  layers. 
This  is  represented  by 


Cl{z)  =  Cl^S{z  -  zi)  +  Cl^S{z  -  Z2)  (3.56) 

Using  this  form  for  the  sifting  property  of  the  delta  function  is  employed  to  remove 
the  integrals  over  z'  and  z".  The  expression  for  the  variance  is  now  reduced  to 


VAR  =  (3.57) 

2.13/:‘‘C^,C^,i"®y  J  du  dv  tnijYnij) 

y-[fi{u,v)  +  2f2{u,v)] 

+2.92k^(7^C^^L~^  J  du  tri(^)/3(u) 

+2.92k^<7‘C^^L~^  J  du  tri(^)/4(ti) 

+(t'‘ 

-  ^lA6k^C-^L~^  J  dwtri(^)/5(u) 

+lA6k~C-^L~^  J  dwtri(^)/6(ti) 

There  are  now  six  functions  (/i  to  /e)  instead  of  four.  This  is  because  the  existence  of 
two  layers  of  turbulence  causes  the  second  and  fourth  integrals  in  equation  3.50  to  each 
split  into  two  integrals,  one  involving  and  one  involving  .  The  six  functions  are  in 
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terms  of  pi  =  1)2  =  Pa  =  Pi  -  Aa:,  and  p^  =  po  -  Ax.  They  are 


fi{u,v)  = 

-2[(p,  +  i)2  +  «2]f  [p2^„2J- 

-2[p?  +  t/2]^  [(pi-i)2  +  r;=]^ 

-2[PI +  «^]^  [(pi  +  T)2  +  t;2j« 

-2  [(pi  -  +  u^]  ^  [pI  +  v^]  ^ 

+  [(P1  +  X)2  +  U^]^  [{P,  +  Lf  +  V^]^ 
^[{p,  +  Lf  +  u^]^  [{p,-Lf  +  v^]^ 
+  [(Pi  -  +  M^]  ‘  [(pi  +  Lf  +  v^]  * 

+  ((Pi  -  Lf  +  vr]  ‘  [(pi  -  Lf  +  v^]  ^ 
+4  [P? +  [pl  +  v^r 
-2[{p,+Lf  +  u^]^  [pI  +  v^-]^ 

-2  (pj  +  U^]  ^  [(P2  -  If  +  v^]  ‘ 

-2  [pj  +  U-]  ^  [(712  +  Lf  +  v^]  ^ 
-2[{p,-Lf  +  n^]^  [pl  +  v^]^ 

+  [(p,+Lf  +  u^]"  [(p,  +  Lf  +  v^‘]^ 
+  [(Pi  +  Lf  +  u^]  ®  [(p,  -  Lf  +  v^]  ® 
+  [(p,-Lf  +  u^]^  [{p,  +  Lf  +  v-]^ 
+  [{p,-Lf  +  u']"  [(p,-Lf  +  v^-y^ 

+4[p;  +  w^]*  [pi+v^]® 

-2[{p2  +  Lf  +  u^]^  [Pl  +  t;=]^ 

-2  [j)l  +  U-]  ‘  [(pi  -  Lf  +  v‘]  ® 

-2  [pI  +  u^]  ^  [(p,  +  Lf  +  V-]  ^ 
-2[{p2-Lf  +  u^]^  [Pl  +  v^-r 
+  [iP2  +  Lf  +  uY  [{Pi  +  Lf  +  v^-]^ 
+  ((P2  +  Lf  +  W-]  ‘  [(Pi  -  Lf  +  ti-]  ^ 
-i-[{P2-Lf  +  u^]^  [(p,  +  Lf  +  v^-]^ 
+  {{P2-Lf  +  u^]"  [(p,-Lf  +  v¥ 


-2  [{p2  +  Lf  +  «-]  ^  [pI  +  v-]  ® 

-2{pI  +  u^]"  l{p,-Lf  +  v^-]i 
-2  [pI  +  ti^] [{p2  +  Lf  +  V-]  * 

-2  ((?;,  -  Lf  +  M-]  ‘  [pI  +  V-]  ® 

+  [{p2  +  Lf  +  ir]  ®  [(yjo  +  Lf  +  v~]  ^ 

+  [(7^2  +  Lf  +  ii"]  ‘  [(7J2  -  Lf  +  v'\  ‘ 

+  [(p2-/.)='  +  tr]"  [(7^3  +  i:)-  +  t;-l^‘ 

+  [(P2  -  Lf  +  W^]  ‘  [(722  -  Lf  +  u*]  ‘ 

The  funccinn  f2{u,v)  has  the  same  form  as  ffu.v)  except  the  pi  is  replaced  bj'  p^  and  tlie 

P2  is  replaced  by  p^.  The  remaining  functions  are 


hin)  = 

2  Irf  +  ^  - !(,),  Kp,  +  /,)>  +  .,=1  > 


(.3.59) 


=  /3(m)  with  721  replaced  by  p2. 


/g(«)  = 

2  (7^3  +  “*]  ^  “  ((7^3  -  T)'  +  «■]  ”  -  ((723  +  Lf  +  «*]  * 


(3.60) 


/elw)  =  U{u)  with  7>3  replaced  by  p^. 


3.5.2  Reduction  of  the  Variance  Expression  In  atmos|)heric  optics,  the  atmospheric 
coherence  diameter,  Jq,  is  a  commonly  encountered  parameter.  In  a  diffraction  limited 
sj'stem  u.sing  a  long  exposure,  the  resolution  will  increase  with  aperture  size.  However, 
the  resolution  will  reach  a  limiting  value  beyond  which  an  incren.se  in  aperture  size  will 
not  affect  icsolution.  The  limiting  value  is  known  as  the  atmospheric  coherence  diameter 
(6:'129-'139).  It  is  often  thought  of  .as  a  measure  of  how  good  the  ‘.seeing"  is  at  any  given 
time.  To  make  the  cxprc.ssion  for  the  v,ari.ance  (and  hence  the  SNR)  more  meaningful,  it 
is  helpful  to  relate  it  to  tq.  Goodmaii  {6:'131)  pro’. ides  an  cxprc.ssion  for  rg-. 


Tg  =  0.1  $.5 


IgQia<l^ 


(3.61) 


Using  the  assumption  of  two-layer  turbulence,  the  integral  expression  in  the  denominator 
may  be  replace  by  Cl^  +  C,;,.  Replacing  A  by  ^  factoring  out  an  L  from  each  of  the 
terms  in  the  integrals,  equation  3.57  may  be  written  in  terms  ol  dimensionless  integrals  as 
follows: 


VAR  = 


(3.62) 


n.98(f)»A' 


£Lr'  r'  f 

cs,\  J-J-, 


(hi'  r/v'  tri(«')tri(i;') 


+6.92(-)^/.-V 

I'll 


+6.92(-)^/,-V 

J’o 


+  J  ^  du  tri(ii')/.;(«') 
du'  tri(«')/s(«') 


r- 

1  +  — 

CL 


[l  +  f  du'  tri(tO/i(«') 

t'o  [  C-^\  J-i 


where  is  the  di.ncnsionless  equivalent  of  /„.  This  form  is  convenient  in  that  and 
now  only  appear  in  a  ra.tio.  It  is  not  necessary  to  assume  exact  vaiii.^s.  only  a  ratio  of 
strengths  at  different  altitudes  in  the  two  layer  model. 


5.5.3  S.^'R  E<iu(tlion  The  SNR  may  now  be  calculated  using  equations  3.47,  3.21, 
and  3.62  in  equation  3.43.  .As  a  final  simplification,  let  /r'  =  aL.  By  using  o'  in  equation 
3.62.  the  factor  L~-  is  common  to  the  numerator  and  denominator  ol  the  SNK  and  is 
thereby  eliminated.  The  equation  for  the  SNR  is 
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da' 


<  3.46 

1 

L 

+3.46 


c- 1"' 

/■’  ,  A 

du  iY\{ll  )f^{u  )  > 

11.98  ^  J  j  du'  dv'  lri(tt')tri(t?') 


+6.92(— 

*'o 


+6.92(— 

I’o' 


+a-(^)- 


1  + 


C? 


T  -1 


'n, 


r’2 

1  +  — 


J  dv  lT\{u')f!,{a') 
J  du'  tri(«')/5(tt') 


3.46 


J  +  ^1  /'  tri(«')/o(»0  ' 


+3.46 


1  + 


c- 

C- 


/•' 

I  ^  f/ti'  lri(H')/i(w') 


In  tiic  ne.xl  cliapler;  .some  lypical  value.s  are  cho.sen  for  tlie  gcomelry  of  figure  3.1  and 
the  SNR  is  c.xamined  a.s  several  parameters  are  allowed  to  vary. 


S.6  Conclasion 

The  correlation  of  wave  front  sloj)e  measurements  has  been  related  to  C'f,  by  an 
integral  e.xpression  involving  a  path  weighting  function.  The  characteristics  of  this  function 
arc  c.xamined  in  the  ne.xl  chapter.  Also,  a  method  to  improve  the  path  weighting  function 
was  discu.ssed  and  an  e.xpression  to  calculate  the  signal-to  noise  ratio  of  the  correlation 
measurement  was  derived.  Thc.sc  are  al.so  c.xamined  in  the  ne.xt  chapter. 


IV.  Findings  and  Results 


4.1  Introduction 

This  chapter  presents  the  findings  and  results  of  this  thesis  research.  The  character¬ 
istics  of  the  unmodified  path  weighting  function  are  discussed  first.  It  is  shown  that  the 
shape  of  the  weighting  function  must  be  modified  to  give  it  the  characteristics  of  a  sampling 
function.  A  method  to  improve  the  shape  of  the  weighting  function  is  presented,  and  the 
resolution  possible  with  the  improved  function  is  calculated.  Finally,  assuming  the  turbu¬ 
lence  is  confined  to  two  layers,  the  signzd-to-noise  ratio  of  the  correlation  measurement  is 
calculated  for  a  sample  case. 


4.2  Characteristics  of  the  Path  Weighting  Function 

In  Chapter  3,  the  path  weighting  function  was  seen  to  be  defined  by  the  equation 


w{z')  = 

-1.46&^£^  j  du'  tri(u')  x 


(4.1) 


if 


(£  +  1)2  +  u'2 


if 


) 


where  z'  is  point  ^llong  the  z  <ixis,  z,  is  the  height  of  the  reference  sources.  Ax  is  the 
separation  of  the  apertures,  Ap  is  the  separation  of  the  reference  sources,  k  =  y,  L  is  the 
dimension  of  the  Lx  L  apertures,  and  p  =  Ap^^  —  Ax  is  the  transverse  separation  of  the 
ray  paths.  Recall  that  if  xu{z')  is  to  be  useful  as  a  sampling  function,  it  should  be  at  a 
maximum  at  the  intersection  point  of  the  ray  paths.  The  function  should  decay  rapidly  to 
zero  away  from  the  intersection  point.  A  function  exhibiting  these  characteristics  could  be 
used  as  a  sampling  function. 
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4.2.1  Unmodified  Path  Weighting  Function  The  normalized  path  weighting  func¬ 
tion  is  plotted  versus  ^  in  figure  4.1.  The  data  points  used  to  create  this  plot  were  generated 
by  the  FORTRAN  computer  program  NORMPATH  contained  in  appendix  A.l. 


Figure  4.1.  Normalized  Path  Weighting  Function 


When  f  =  0,  the  ray  paths  are  at  a  point  of  intersection  and  the  function  is  at  a 
maximum.  As  ^  increases,  the  ray  paths  diverge  and  the  function  falls  off.  Since  the 
expression  for  the  path  weighting  function  was  made  dimensionless,  the  curve  in  figure  4.1 
is  the  general  shape  for  all  values  of  L.  For  example,  ML  =  Ivi,  the  values  on  the  horizontal 
axis  range  from  0  to  100m.  If  i,  =  1cm,  the  values  range  from  0  to  Im.  A  change  in  the 
value  of  L  is  therefore  seen  to  control  the  decay  time  of  the  curve.  As  evidenced  by  figure 
4.1,  the  path  weighting  function  does  exhibit  the  general  required  sampling  charac. eristics. 
It  has  a  maximum  value  at  ^  =  0  and  falls  off  as  ^  increases.  However,  the  fiiiictio.i  docs 
not  fall  off  rapidly  and  never  reaches  a  value  of  zero.  The  most  rapid  decay  occurs  when 
L  is  smallest,  but  the  function  still  does  not  approach  zero  quickly  enough.  Therefore, 
the  path  weighting  function  must  be  modified  to  more  closely  approximate  a  sampling 
function,  as  shown  in  the  next  section. 
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4.2.2  Modified  Path  Weighting  Function  Tiic  method  used  to  improve  the  shape  oj 
the  path  weighting  function  is  to  combine  weighting  functions  generated  by  two  apertures 
of  sizes  Li  and  L-^-  In  particular,  this  should  remove  the  common  component  of  the 
average  tilt  and  yield  a  function  with  a  sharp  falloff.  By  subtracting  a  weighting  function 
generated  by  an  aperture  of  size  L2  from  that  generated  by  an  aperture  of  size  Li.  the 
common  low  spatial  frequency  components  arc  removed.  These  low  frequency  components 
have  a  large  correlation  over  large  separations  of  the  ray  paths.  The  large  correlation  values 
associated  with  these  low  frequency  components  are  what  caused  the  previous  results  to 
be  unacceptable.  In  effect,  the  subtraction  is  equivalent  to  a  high  pass  filtering  of  the 
weighting  function.  A  plot  of  the  modified  path  weighting  function  is  contained  in  figure 
4.2  for  three  ratios  of  The  data  points  used  to  create  this  plot  were  generated  bv-  the 
FOIITR.A.N  computer  program  MODPWF  contained  in  appendix  A.2. 


Figure  4.2.  Modified  Path  W^ghling  Function 


As  with  the  normalized  path  weighting  function  plot,  this  method  yields  a  family  of 
curves.  The  horizontal  spread  is  controlled  by  the  size  of  Li  and  Lo-  As  Li  and  Xo  get 
larger,  the  .spread  increases.  This  is  because  the  low  frequency  components  are  correlated 
over  a  larger  area.  This  action  introduces  a  larger  common  tilt  component.  The  rate  of 
fallofF  to  zero  is  controlled  by  the  ratio  As  the  ratio  approaches  1,  decays  more 

rapidly.  Therefore,  the  best  sampling  characteristics  are  exhibited  by  a  curve  that  has 
small  values  for  Li  and  latio  that  approaches  1.  If  Li  =  lc77t  and  Xj  =  2cvi, 

the  horizontal  axis  in  figure  4.2  is  scaled  from  0  to  Im.  For  a  ratio  of  0.9,  it  is  seen  that 
the  path  weighting  function  goes  to  0  when  the  transverse  ray  path  separation  is  only  0.5 
meters,  which  is  very  close  to  the  intersection  point.  The  plot  shows  that  the  modified 
path  weighting  function  exhibits  excellent  sampling  characteristics.  The  function  falls  off 
rapidly  and  stays  at  a  zero  value  and  could  therefore  be  used  to  extract  the  value  of  C;^ 
from  the  correlation  value. 

4.3  Rest  ion  of  (he  Ruth  Weighling  Function 

It  has  been  shown  that  the  modified  path  weighting  function  can  be  used  as  a  sam¬ 
pling  function.  It  now  remains  to  determine  the  resolution  of  the  function  as  it  relates 
to  the  measurement  geometry.  The  width  of  the  function  w(^')  may  be  interpreted  as  a 
measure  of  the  vortical  resolution.  The  narrower  iv{z')  is,  the  more  accurate  will  be  the 
value  of  that  is  calculated  from  C,(A.r).  This  may  be  directly  related  to  the  sifting 
property  of  the  Dirac  della  function. 

To  measure  the  vertical  resolution,  the  method  used  by  Welsh  (18)  is  adopted.  Recall 

that 

z' 

— A.t  (4.2) 


by 


At  the  intersection  of  the  ray  paths,  p  =  0  and  the  inler.scction  altitude  is  given 


Let  a  measure  of  the  width  of  tu(z')  be  given  by  the  point  p  =■  p'  where,  for  example, 
p'  is  the  e“*  point  of  w(^).  Designate  the  altitude  corresponding  to  this  point  as  Zp>. 
Solving  for  Zp-: 


(p'  +  Ax)  z. 


(4.4) 


The  width  of  w{z),  with  respect  to  z^  is 


^Z  =  Zq-  Zp 


P'  z. 


(4.5) 


Let  0  designate  the  angle  between  ray  paths.  As  seen  from  figure  3.1,  this  is  also  the 
angular  separation  of  the  point  sources.  0  may  be  expressed  as 


^=arctan(^)  (4.6) 

Solving  for  Apa,, 


Ap^  =  22:,tan(^)  (4.7) 


Using  equation  4.7,  equation  4.5  may  be  expressed  cis 


A^  = 


P' 

2tan(|) 


(4.8) 
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As  an  example,  referring  to  figure  4.2,  for  Li  =  1cm,  L2  =  2cm,  and  ^  =  0.5,  it  is 
seen  that  the  width  of  is  approximately  =  0.8m.  Using  this  value  for  p',  vertical 
resolution  is  plotted  as  a  function  of  9  in  figure  4.3. 


Figure  4.3.  Vertical  Resolution 


The  plot  clearly  shows  that  resolution  is  related  to  separation  of  the  point  sources. 
For  example,  if  10  meter  resolution  is  required,  9  must  be  approximately  0.4  degrees.  In 
general,  resolution  increases  as  source  separation  increcises. 
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/iJi  Signal-lo-Noise  Ratio  of  the  Correlation  Measurement 

In  order  to  quantify  the  usefulness  of  the  methods  presented  in  this  thesis,  the  SNR 
of  the  correlation  measurement  is  calculated.  The  reader  is  referred  to  Chajrter  .3  for  the 
derivation  and  defining  equations.  It  is  important  to  note  the  SNR  analysis  is  valid  for  the 
case  of  the  unmodified  path  weighting  function  only.  While  this  provides  a  good  initial 
look  at  the  problem,  the  analysis  should  be  e.xtended  in  the  future. 

To  evaluate  the  SNR  for  a  typical  case,  it  was  necessary  to  choose  values  for  several 
parameters.  The  reader  will  recall  that  a'~  x  was  a  parameter  in  the  SNR  equation 

that  characterizes  the  noise  in  the  process.  This  parameter  will  be  allowed  to  vary  to  see 
how  noise  effects  the  SNR.  The  structure  constant  ratio  ^  =  2.5  was  used  based  on  the 

"I 

data  in  figure  2.1.  Assuming  the  turbulence  is  confined  to  two  layers,  the  peaks  of  the 
curve  at  1  and  10/cm  are  chosen.  The  approac.h  \.ken  is  to  divide  the  area  under  the  curve 
into  two  sections  and  assume  it  exists  as  delta  functions  at  the  peak  values.  The  area 
under  the  curve  from  1  to  10km  is  calculated  and  assumed  to  exist  as  a  delta  function  at 
Ikm.  In  the  same  fashion,  the  area,  from  10  to  20/cm  i.s  calculated  and  assumed  to  exist  as 
a  della  function  at  20/cm.  The  resulting  ratio  is  2.5. 

The  reference  sources  are  assumed  to  be  laser  guidestars  at  an  altitude  of  c,  =  100km. 
Guidestar  separation,  A/;,,  is  1/cm.  The  data  points  used  to  create  the  plots  in  the  following 
subsections  were  generated  by  the  FORTRAN  computer  program  S.NR.  in  appendix  A. 3. 


/i.jj.l  SNR  as  a  Function  of  Noise  Welsh  and  Gardner  (19:1919)  provide  an  ex¬ 
pression  for  the  tilt  measurement  error  (noise)  as  follows: 


(T  = 


0.867777 

iVar-Q 


L  >  7-0 


(4.9) 


where  77  is  a  parameter  used  to  account  for  imperfections  in  the  detector  array  and  N  is 
the  total  subaperture  photon  count.  77  =  1.5  is  used  as  a  typical  value. 

Using  equation  4.9  and  recalling  o'  —  aL,  it  is  possible  to  relate  a'-  to  the  photon 
count  for  the  values  used  in  the  sample  case.  The  resulting  equation  is 


N  _  (0.867777)2 


(4.10) 


where  the  parameter  has  been  introduced  to  coincide  with  the  terms  in  the  SNR. 

equation  in  Chapter  3.  The  photon  count  may  now  be  related  to  the  SNR  using  equation 
4.10.  This  relationship  is  plotted  in  figure  4.4. 
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The  plot  shows  tliat  a  photon  count  of  greater  llian  200  produces  an  SNR  very  close 
to  the  noiseless  value.  Assuming  that  bright  guidestars  are  used  (i.c.  created  with  sufficient 
power),  photon  noise  should  not  be  a  concern.  The  remainder  of  the  SNR  analysis  will 
therefore  assume  a  noiseless  system  (i.e.  a'  =  0). 

SNR  as  a  Function  of  Ray  Path  Separation  In  this  subsection,  the  effect  on 
the  SNR  of  changing  the  ray  path  separation  is  investigated.  There  are  two  ways  to  change 
the  path  separation;  by  changing  the  separation  of  the  subapertures  or  by  changing  the 
separation  of  the  point  sources.  In  figure  d.5,  the  subaperture  separation,  Aa:,  is  allowed 
to  vary.  Recall  that  Zi  is  the  height  of  the  first  turbulence  layer,  is  the  height  of  the 
second  turbulence  layer,  z,  is  the  guidestar  height,  and  Apr  is  the  guidestar  separation. 
The  separation  of  ray  paths  is  given  by  Apr^  -  Ax.  Therefore,  the  two  points  of  interest 
on  the  plot  are  when  A.(;  =  Ap^^  and  A.’C  =  Ap^^.  These  points  correspond  to  a  path 
intersection  point  for  each  turbulence  layer  and  occur  at  Ax  =  10  and  Ax  =  100.  In  figure 
4.6,  the  guidestar  separation,  Ap^.  is  allowed  to  vary.  The  corresponding  intersection  points 
for  the  turbulence  layers  are  at  Ap^  =  lOOrrj  and  Ap*  =  Ihm. 


I-!) 
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Figure  4.6.  SNR  as  a  Function  of  Source  Separation 
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As  seen  from  the  plots,  the  SNR  reaches  two  maxiimiin  values.  While  the  plots  ap])ear 
to  be  in  reverse  order,  it  is  important  to  realize  that  the  SNR  values  are  consistent  between 
the  two  plots.  For  e.\ainple,  the  value  of  0.259  at  A.r  =  lOm  and  =  14:7/1  is  consistent 
between  the  two  plots.  This  becomes  clear  when  one  considers  what  is  happening  when 
either  the  source  or  subaperture  separation  is  increased.  The  net  effect  is  to  change  ray 
path  separation,  however,  the  change  is  initiated  at  opposite  ends  and  this  causes  the  plots 
to  appear  reversed. 

The  difference  in  the  height  of  the  peaks  in  the  turbulence  layers  may  be  directly 
related  to  the  strength  of  the  turbulence  at  the  layers.  The  stronger  turbulence  is  causing 
a  higher  correlation  value  and  hence  a  larger  SNR.  It  may  be  assumed  that  it  will  be  easier 
to  detect  stronger  turbulence  values  with  this  method. 

The  variation  in  the  plot  of  the  SNR  indicates  that  different  combinations  of  .source 
and  subaperture  separations  yield  different  values  of  the  SNR.  Depending  on  the  applica¬ 
tion  and  physical  limitations  of  an  e.\periment  (i.e.  physical  space  available  to  change  the 
separation  of  the  wave  front  sensors),  one  may  adjust  the  interrelationship  for  an  optimum 
SNR. 


4.4.3  Improvement  of  the  SNR  The  SNR  plots  in  the  previous  section  did  not 
ex’ceed  a  value  of  1.  For  the  method  to  be  useful,  it  is  necessary  to  consider  how  the  SNR 
may  be  increa.scd.  'Phe  reader  wdll  recall  an  initial  assuin|)tion  in  Chapter  3  restricted 
the  wave  front  sensors  to  a  single  subaperture.  If  sensors  vvith  multiple  subaperliires  are 
used,  an  array  of  (orrelation  values  will  be  obtained.  This  is  done  by  pairing  all  po.ssiblc 
combinations  of  the  subapertures,  each  pair  being  separated  by  A.r.  U.sing  an  array  of  n 
values  should  improve  the  SNR  by  a  facior  of  n.  This  action  should  raise  the  SNR  to  a 
level  suitable  for  detection. 

4. -5  Ohlciritu/  a  Vertical  Rrojilc 

As  was  .seen  in  the  previous  section,  changing  the  subaperture  .separation  moves  the 
path  inter.scction  point  vertically.  This  feature  may  be  e.\i)loiled  to  obtain  a  \ertical  profile 
ofC’,^.  If  sen.sors  with  multiple  subapertures  are  u.sed.  it  is  pos.sible  u.se  many  different  values 
of  A.7:.  Fairing  all  po.s.sible  combinations  of  subapei  t  iircs  separated  by  all  values  of  A.r  will 
result  in  an  array  of  correlation  values  that  correspond  to  different  intei.section  altitudes. 
The.se  values  may  be  exploited  to  extract  the  value  of  C‘  at  the  different  altitudes. 


4- 6  Conclusion 

It  has  been  shown  that,  using  the  subtraction  method,  the  path  weighting  function 
can  be  modified  to  approximate  a  sampling  function.  C,  can  be  directly  related  to  by 
the  modified  function.  The  resolution  plot  indicates  that  source  separation  determines  the 
vertical  resolution.  As  source  separation  increases,  vertical  resolution  increeises.  The  SNR 
calculations  for  the  sample  case  show  that  wave  front  sensors  with  multiple  subapertures 
will  be  necessary.  If  multiple  values  are  obtained  to  raise  the  SNR,  it  will  be  possible  to 
extract  the  value  of  from  the  value  of  C,.  Using  multiple  subapertures  also  allows  one 
to  obtain  a  vertical  profile. 
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V.  Conchsions  and  Recommendations 


5.1  Overview 

This  chapter  contains  the  conclusions  drawn  from  this  research.  It  also  contains  the 
authors  recommendations  for  further  study. 

5.2  Conclusions 

This  research  has  shown  that  the  spatial  correlation  of  wave  front  slopes  m<ay  be  used 
to  calculate  the  value  of  the  atmospheric  structure  constant.  The  correlation  process  yields 
an  integral  exp.''' ^.-jion  relating  C"’  to  the  correlation  value  by  means  of  a  path  weighting 
function.  It  is  possible  to  shape  the  weighting  function  into  a  sampling  function  and  to 
use  the  sifting  property  of  the  sampling  function  to  directly  relate  CJ,  to  the  correlation 
measurement. 

It  is  difficult  to  compare  the  accuracy  of  this  method  to  other  methods.  Typically, 
other  methods  of  calculating  Cl  arc  experimentally  implemented  and  the  results  of  the 
experiment  are  compared  against  known  values  of  Cl.  That  has  not  been  done  yet  with 
this  research.  However,  based  on  the  sample  SNR  calculation  for  two-layer  turbulence,  it  is 
anticipated  this  method  will  yield  accurate  results.  Full-scale  modeling  «nd  testing  sitould 
prove  this  assertion  to  i)e  true. 

5.3  Recommendations 

The  following  suggestions  for  further  study  are  made. 

1.  The  S.NR  analysis  should  be  extended  beyond  the  a.ssumption  that  ihe  turbulence 
is  conrined  to  two  layers.  Increasing  the  number  of  turbulent  lavers  will  more  closel\ 
approximate  a  real-life  situation. 

2.  'I'he  SNR.  anal\.sis  should  be  done  for  the  case  of  the  modified  path  weighting  function. 
.3.  Using  an  accepted  model  for  atmospheric  turbulence.,  this  re.search  method  should  be 
computer  modeled.  The  spatial  correlation  of  \\a\e  front  phase  slopes  should  be  calculated 
for  the  modeled  scenario  and  compared  to  known  lesult.s  to  determine  the  accuraev  of  this 
method. 

4.  This  method  should  be  experimentalise  verified  and  the  results  should  be  compared  to 
those  obtained  bv  other  mea.surement  methods. 


5.Jt  Summary 

The  spatial  correlation  of  wave  front  slopes  from  two  point  sources  theoretically  yields 
a  value  for  C^.  Experimental  verification  of  the  method  will  allow  it  to  become  another 
tool  to  calculate  the  magnitude  of  atmospheric  turbulence. 
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Appendix  A.  Compuler  Programs 


This  appendix  contains  the  FORTRAN  computer  code  used  to  obtain  tlie  plots  in 
Chapter  4  of  this  thesis.  The  code  wa.s  written  by  the  author  with  the  exception  of  QSIMP, 
TRAPZD,  and  INTEG2D.  These  were  adapted  from  the  book  Numerical  Recipes:  The  Art 
of  Scientific  Computing  by  William  H.  Press.  Brian  P.  Flannery,  Saul  'V.  Teukolsky,  and 
William  T.  Velterling  published  by  Cambridge  University  Press,  New  York,  1986,  jiages 
110-1.30. 


AJ  Program  NORM PATH 


PROGRAM  NORMPATH 

* 

*  This  program  computes  the  values  for  the 

*  normalized  path  weighting  function. 

* 

EXTERNAL  KTFUNC 
COMMON  /PARAMS/P,L 
REAL  P,PPRIHE,L, NORM, STEP 
INTEGER  I 

100  FORMAT  (F10.3) 

OPEN(UNIT=10,  NAME='PR0G1.DAT’ ,  STATUS= ’ NEW » ) 

* 

READ(5,100)  L 
READ (5, 100)  STEP 
P=0. 

* 

CALL  QSIMP(WTFUNC,-1,1,S) 

N0RM=S*(L**(-l./3.)) 

DO  110  1=0,100 
PPRIME=P/L 

CALL  QSIMP (WTFUNC,-1,1,S) 
S=(S/N0RM)*(L**(-l./3.)) 

WRITE(10,100)  PPRIME.S 
P=P+STEP 
no  CONTINUE 
* 

CLOSE(UNIT=10) 

STOP 

END 

^  ^  4:  :4c  sfr  4c  4c  4c  4c  ;|c  4c 

REAL  FUNCTION  WTFUNC(X) 
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COMMON  /PARAMS/P,L 
REAL  P,L 

WTFUNC=-Cl.-ABS(X))*( 

C  2.*(((P/L)**2.+(X)**2.)**(5./6.))- 

C  ( ( (P/L) - 1 . ) **2 . + (X) **2 . ) ** (5 . /6 . ) - 

C  (((P/L)+l.)**2.+(X)**2.)**(5./6.) 

C  ) 

RETURN 

END 

SUBROUTINE  TRAPZD(FUNC,A.B,S,N) 

IF  (N.EQ.l)  THEN 
S=0 . 5* (B-A) * (FUNC ( A) +FUNC(B) ) 

IT=1 

ELSE 

TNM=IT 

DEL=(B-A)/TNM 

X=A+0.5*DEL 

SUM=0. 

DO  200  J=1,IT 
SUM=SUM+FUNC(X) 

X=X+DEL 
200  CONTINUE 

S=0 . 6* (S+ (B-A) *SUM/TNM) 

IT*2*IT 

ENDIF 

RETURN 

END 

★  ♦***!l<*******>l'**************j(t**********j|£****!it***** 

SUBROUTINE  QSIMP(FUNC,A,B,S) 

PARAMETER  (EPS=l.E-5,  JMAX=20) 

0ST=-i.E30 
0S=  -1.E30 
DO  300  J=1,JMAX 

CALL  TRAPZD(FUNC,A,B,ST,J) 

S=(4.*ST-0ST)/3. 

IF  (ABS(S-OS).LT.EPS*ABS(OS))  RETURN 
OS=S 
OST=ST 
300  CONTINUE 

PAUSE  ’Too  many  steps.’ 

END 


A.2  Program  MODPWF 


PROGRAM  MODPWF 

* 

*  This  progrean  subtracts  a  weighting  function 

*  for  an  aperture  of  size  L2  from  a  weighting 

*  function  for  an  aperture  of  size  Li. 

*  The  values  are  normalized  for  a  maximum  of  1. 

* 

EXTERNAL  HTFUNCL1,HTFUNCL2 
COMMON  /PARAMS/P,L1,L2, RATIO 
REAL  P,PPRIME,L.L1.L2,RATI0,N0RM,STEP 
INTEGER  I 

100  FORMAT  (F10.3) 

OPEN(UNIT=10,  NAME=' PR0G1.DAT',  STATUS='NEW») 

* 

P-=0. 

READ (5, 100)  LI 
READ (5, 100)  RATIO 
READ (5, 100)  STEP 
L2=L1/RATI0 

* 

DO  110  1=0,100 
L=L1 

CALL  QSIHP(MTFUNCL1,-1.,1.,SA) 
SA=SA*ai**(-l./3.)) 

L=L2 

CALL  QSIMPCWi. JNCL2,-1.,1.,SB) 
SB=SB*(L2**(-l./3.)) 

S=SA-SB 

IF  (I  .EQ.  0)  THEN 
H0RM=S 
EHDIF 
S*S/N0RM 
PPRIME=P/L1 

WRITEClO.lOO)  PPRIME.S 
P=P+STEP 
110  CONTINUE 
♦ 

CL0SE(UNIT=10) 

STOP 

END 

^i**^^^c*iti*,^**t^i********^i^**:HL*t**********t*********^* 

REAL  FUNCTION  WTFUHCLl(X) 

COMMON  /?ARAMS/P,L1,L2, RATIO 


REAL  P, LI, L2, RATIO 
WTFUNCL1=- ( 1 . -ABS (X) ) * ( 

C  2.*(((P/Ll)**2.+(X)**2.)**(5./6.))- 

C  ((P/Ll-l)**2.+(X)**2.)**(5./6.)- 

C  ((P/Li+l)**2.  +  (X)**2.)=i”i'(5./6.) 

C  ) 

RETURN 

END 

REAL  FUNCTION  WTFUNCL2(X) 

COMMON  /PARAMS/P,L1,L2, RATIO 
REAL  P, LI, L2, RATIO 
WTFUNCL2=- ( 1 . -ABS (X) ) * ( 

C  2.*(((P/L2)**2.  +  (X)*>i<2.)**(5./6.))- 

C  ((P/L2-1)**2.+(X)**2.)**(5./6.)- 

C  ((P/L2+1)**2.+(X)**2.)**(5./6.) 

C  ) 

RETURN 

END 

SUBROUTINE  TRAPZD(FUNC,A,B,S,N) 

IF  (N.EQ.l)  THEN 
S=0 . 5* (B-A) * (FUNC (A) +FUNC(B) ) 

IT=1 

ELSE 

TNM*IT 

DEL=(B-A)/TNM 

X=A+0.5i<DEL 

SUM=0. 

DO  200  J=1,IT 
SUM=SUM+FUNC(X) 

X=X+DEL 
200  CONTINUE 

S=0 . 5* (S+ (B-A) *SUM/TNM) 

IT=2*IT 
ENDIF 
'  RETURN 
END 

SUBROUTINE  QSIMP<FUNC,A,B,S) 

PARAMETER  (EPS=l.E-5,  JMAX=20) 

0ST=-l.E3O 
0S=  -1.E30 
DO  300  J=1,JMAX 

CALL  TRAPZD(PUNC,A,B,ST,J) 
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S=(4.*ST-0ST)/3. 

IF  (ABS(S-OS).LT.EPS*ABS(OS))  RETURN 
OS=S 
OST=ST 
300  CONTINUE 

PAUSE  ’Too  many  steps.’ 

END 
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A. 8  Program  SNR 


PROGRAM  SNR 

* 

*  This  program  computes  the  mean  and  variance 

*  of  the  correlation  of  wave  front  phase  slopes. 

*  It  then  calculates  the  signal-to-noise  ratio 

*  from  the  relationship  SNR=MEAN/SQRT(VAR) . 

* 

EXTERNAL  WTFUNCAl ,WTFUNCA2 .WTFUNCB 
EXTERNAL  WTFUNCC , WTFUNCD , WTFUNCE 
EXTERNAL  G,H 

COMMON  /PARAMS/P1,P2,P3,P4 

COMMON/FLGS/FLAG 

REAL  P1,P2,P3,P4 

REAL  SA,SB,SC,SD,SE,S 

REAL  CA,CB,CC,CD,CE 

REAL  L_R0,CN21 .DELTAX ,DEL_PX ,N0IVAR 

REAL  MEAN, VAR, SNR 

INTEGER  I, FLAG 

* 

100  FORMAT  (F10.3) 

0PEN(UNIT=10,  NAHE=’ PR0G1.DAT’,  STATUS= ' NEW ’ ) 

* 

L_R0=10, 

CN21=2.5 

DELTAX=10. 

DEL_PX=1000. 

P1=DEL_PX/100. 

P2=DEL_PX/10. 

P3=(P1-DELTAX) 

P4=(P2-DELTAX) 

N0IVAR=O. 

* 

DO  no  1=0,100 
CA=11.98*(1/(CN21+1/CN21)) 
CB=6.92*(L_R0**(-5./3.))* 

C  (1/(1+CN21))*N0IVAR 

CC=6.92*(L_R0**(-5./3.))* 

C  (1/(1+1/CN21))*M0IVAR 

CD=-3.46*(1/(1+CN21)) 
CE=-3.46*(1/(1+1/CN21)) 

FLAG=1 

CALL  INTEG2D(-1. ,1. ,S) 

SA1=S 


.A-C) 


FLAG=2 

CALL  INTEG2D(-1. ,1. ,S) 

SA2=S 

CALL  QSIMP(WTFUNCB,-1..1-.SB) 

CALL  QSIMP(WTFUNCC,-1.,1.,SC) 

CALL  QSIMP(WTFUNCD,-l.,i.,SD) 

CALL  QSIMP(WTFUNCE,-1..1.,SE) 

MEAN=ABS (CD*SD+CE*SE) 

VAR= (CA* (SA1+SA2) ) + (CB*SB) + (CC*SC) + 

C  (N0IVAR**2 . *L_RO** (- 10 . /3 . ) - (MEAN**2 , ) 

SNR=MEAN/SQRT(VAR) 

WRITE (10, 100)  SNR 

* 

*  Depending  on  the  parameter  being  modified  to  see 

*  the  effect  on  the  SNR,  code  must  be  added  here. 

*  For  example,  if  the  effect  of  noise  on  the  SNR  is 

*  the  quantity  of  interest,  N0IVAR=N0IVAR+1.  may 

*  be  used. 

* 

110  CONTINUE 
* 

CL0SE(UNIT=10) 

STOP 

END 

SUBROUTINE  INTEG2D(X1,X2,S) 

EXTERNAL  H 

COMMON  /PARAMS/P1,P2,P3,P4 

REAL  P1,P2,P3,P4 

CALL  qSIMPX(H,-l.,l.,S) 

RETURN 

END 

*:)t:t!**!t!**J(!J(t**>(>****lt<*l(<*Jt<*****J|'********J|t**!|C***;(C!t!****** 

REAL  FUNCTION  H(XX) 

EXTERNAL  G 

C0MM0N/XY/X,Y 

COMMON  /PARAMS/P1,P2,P3,P4 

REAL  P1,P2,P3,P4 

X=XX 

CALL  QSIMPY(G,-1.,1.,S) 

H=S 

RETURN 

END 

REAL  FUNCTION  G(YY) 
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EXTERNAL  WTFUNCAl ,WTFUNCA2 

COMMON/XY/X,Y 

C0MH0N/PARAMS/P1,P2,P3,P4 

COMMON/FLGS/FLAG 

INTEGER  FLAG 

REAL  P1,P2,P3,P4 

Y=YY 

IF  (FLAG  .EQ.  1)  THEN 

G=WTFUNCA1(X,Y) 

ELSE 

G=WTFUNCA2(X,Y) 

ENDIF 

RETURN 

END 

^  sfe  ^  %  )|e  Jtcfc  }|(  Jt:  ^  ^  3k  ¥  ^  ^  ^  ^  3f(  3|c  :|c  ^  ^  ^  ^  ^  % 

REAL  FUNCTION  WTFUNCAl (X.Y) 

COMMON  /PARAMS/P1,P2,P3,P4 
REAL  P1.P2,P3,P4 

WTFUNCA1= ( 1 . - ABS (X) ) * ( 1 . -ABS (Y) ) * ( 

C  4.*(Pl**2.+X**2.)**(S./6.)* 

C  (Pl**2.+Y**2.)**(5,/6.)+ 

C  ((Pl+l.)**2.+X**2.)**(5./6.)* 

C  ((Pl+l.)**2.+Y**2.)**(S./6.)+ 

C  ((Pl-l.)**2.+X**2.)**(5./6.)* 

C-  ((Pl+l.)**2.+Y**2.)**(5./6.)+ 

C  ((Pl+l.)**2.+X**2.)**(5./6.)* 

C  ((Pl-l.)**2.+Y**2.)**(5./6.)+ 

C  ((Pl-l.)**2.+X**2.)**(5./6.)* 

C  (<Pl-l.)**2.+Y**2.)**(5./6.)- 
C  2.*((Pl+l.)**2.+X**2.)**(5./6.)* 

C  ((Pl)**2.+Y**2.)**(5./6.)- 
C  2.*((Pl-l.)**2.+X**2.)**(5./6.)* 

C  ((Pl)**2.+Y**2,)**(5./6.)- 
C  2.*((P1)**2.+X**2,)**(5./6.)* 

C  ((Pl-l.)**2.+Y**2.)**(5./6.)- 
C  2.*((Pl)**2,+X**2.)**(5./6.)* 

C  ((Pl+l.)*>i'2.+Y**2.)**(5./6.)+ 

C  4.*(Pl**2.+X*>i‘2.)**(5./6.)* 

C  (P2>t<*2.+Y**2.  )**(5./6.)  + 

C  ((Pl+l.)**2.+X**2.)**(5./6.)* 

C  ((P2+l.)**2.+Y**2.)**(5./6.)+ 

C  ((Pl-l.)**2.+X**2.)**(5./6.)* 

C  ((P2+l.)**2.+Y**2.)**(5./6.)+ 

C  ((Pl+l.)**2.+X**2.)**(5./6,)+ 

C  (CP2-l.)**2.+Y**2.)**(5,/6.)+ 
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C  ((Pl-l.)**2.+x>t<*2.)**(5./6.)>t< 

C  ((P2-l.)**2.+Y**2.)**(5./6.)- 
C  2.*((Pl+l.)*>i‘2.+X**2.)**(5./6.)* 
C  ((P2)**2.+Y**2.)**(5,/6.)- 
C  2.*((Pl-l.)**2.+X**2.)**(5./6.)* 
C  ((P2)**2.+Y**2.)**(5,/6.)- 
C  2.*((Pl)**2.+X**2.)**(5./6.)* 

C  C(P2-l.)**2.+Y**2.)**(5./6.)- 
C  2.*((Pl)**2.+X**2.)**(5,/6.)* 

C  ((P2+l.)**2.+Y**2.)**(5./6.)+ 

C  4.*(P2**2.+X**2.)**(5./6.)* 

C  (Pl**2.+Y**2.)**(5./6.)+ 

C  C(P2+l.)**2.+X**2.)**(5./6.)* 

C  ((Pl+l.)**2.+Y**2.)**(5./6.)+ 

C  ((P2-l.)**2.+X**2.)**(5./6.)* 

C  ((Pl+l.)**2.+Y**2.)**(5./6.)+ 

C  ((P2+l.)**2.+X**2.)**';5./6.)* 

C  ((Pl-l.)**2.+Y**2.)*-i'(5./6.)  + 

C  ((P2-l.)**2.+X**2.)**(5./6.)* 

C  ((Pl-l.)**2.+Y**2.)**(5./6.)- 
C  2 . * ( CP2+1 . ) **2 . +X**2 . ) ** (5 . /6 . ) * 

C  C(Pl)**2.+Y**2.)**(5./6.)- 
C  2.*((P2-1.)*>k2.+X**2.)**(5./6.)* 
C  ((Pl)**2.+Y**2.)**C5./6.)- 
C  2.*(CP2)>k*2.+X**2.)**(5./6.)* 

C  ((Pl-l.)**2.+Y**2.)**(5./6.)- 
C  2.*((P2)**2.+X**2.)**(5./6.)* 

C  ((Pl+l.)**2.+Y**2.)**(5./6.)+ 

C  4.*(P2**2.+X**2.)**(5./6.)* 

C  (P2**2.+Y**2.)**C5./6.)+ 

C  ((P2+l.)**2.+X**2.)**(5./6.)* 

C  ((P2+l.)**2.+Y**2.)**(5./6.)+ 

C  ((P2-l.)**2.+X**2.)**(5./6.)* 

C  ((P2+l.)>K*2.+Y**2.)**(5,/6.)  + 

C  ((P2+l.)**2.+X**2.)**(5./6.)* 

C  ((P2-l.)**2.+Y**2.)**(5./6.)+ 

C  ((P2-l.)**2.+X**2.)**(5./6.)* 

C  ((P2-l.)**2.+Y**2.)**(5./6.)- 
C  2.*((P2+l.)**2.+X**2.)**(5./6.)* 
C  ((P2)**2.+Y**2.)**(5,/6.)- 
C  2.*((P2-l.)**2.+X**2,)**(5./6.)* 
C  ((P2)**2.+Y**2.)**(5./6.)- 
C  2.*((P2)**2.+X**2.)**(5./6.)* 

C  ((P2-1.)**2.+Y**2.)**(5./6.)- 
C  2.*((P2)**2.+X**2.)**(5,/6.)* 


A-9 


C  ((P2+l.)**2.+Y**2.)**(5./6.) 

C  ) 

RETURN 

END 

REAL  FUNCTION  WTFUNCA2(X,Y) 

COMMON  /PARAMS/P1,P2,P3,P4 

REAL  P1,P2,P3,P4 

WTFUNCA2= ( 1 . -ABS (X) ) * ( 1 . -ABS (Y) ) *2* ( 
C  4.*(P3**2.+X**2.)**(5./6.)* 

C  (P3**2.+Y**2.)**(5./6.)+ 

C  ((P3+l.)**2.+X**2.)**(5,/6.)* 

C  ((P3+l.)**2.+Y**2.)**(5./6.)+ 

C  ((P3-l.)**2.+X**2.)**(5./6.)* 

C  ((P3+l.)**2,+Y**2.)**(5./6,)+ 

C  ((P3+l.)**2.+X**2.)**(5./6.)* 

C  ((P3-l.)**2.+Y**2.)**(5./6.)+ 

C  ((P3-l.)**2.+X**2.)**(5./6.)* 

C  ((P3-l.)**2,+Y**2.)**(5./6.)- 
C  2.*((P3+l.)**2.+X*>t‘2.)**(5./6.)* 
C  ((P3)**2.+Y**2.)**(5./6.)- 
C  2.*((P3-l.)**2.+X**2.)**(5,/6.)* 
C  ((P3)**2.+Y**2.)**(5./6.)- 
C  2.*((P3)**2.+X**2.)**(5./6.)* 

C  ((P3-l.)**2.+Y**2.)**(5./6.)- 
C  2.*((P3)**2,+X**2.)**(5./6.)* 

C  ((P3+l.)**2.+Y**2.)**(5./6.)+ 

C  4.*(P3**2.+X**2.)**(5./6.)* 

C  (P4**2.+Y**2.)**(5./6.)+ 

C  ((P3+l.)**2.+X**2.)**(5./6.)* 

C  <(P4+l.)**2.+Y**2.)**(5./6.)+ 

C  ((P3-1.)**2.+X**2.)**C5./6.)* 

C  ((P4+l.)**2.+Y**2.)**(5./6.)+ 

C  ((P3+l.)**2.+X**2.)**(5./6.)* 

C  ((P4-l.)**2.+Y**2.)**(5./6.)+ 

C  ((P3-1.)**2.+X**2.)**(5./6.)* 

C  ((P4-l.)**2.+Y**2.)**(5./6.)- 
C  2.*((P3+l.)**2.+X**2.)**(5./6.)* 
C  ((P4)**2,+Y**2.)i‘*(5./6.)- 
C  2.*((P3-l.)**2.+X**2.)**(5./6.)* 
C  ((P4)**2.+Y**2.)**(5./6.)- 
C  2.*((P3)**2.+X**2.)**(5./6.)* 

C  ((P4-1.)**2.+Y*>k2.)**(5./6.)- 
C  2,*((P3)>k*2,+X**2.)**(5./6.)* 

C  ((P4+l.)**2.+Y**2.)**(5./6.)+ 
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C  4.*(P4>t<*2.+X**2.)**(5./6.)* 

C  (P3**2.+Y**2.)**(5./6.)+ 

C  ((P4+l.)**2.+X**2.)**(5./6.)* 

C  ((P3+l.)**2.+Y**2.)**(5./6.)+ 

C  ((P4-l.)**2.+X**2.)**(5./6.)* 

C  ((P3+l.)**2.+Y**2.)**(5./6.)+ 

C  ((P4+l.)**2.+X**2.)**(5./6.)* 

C  ((P3-l.)>t'*2.+Y**2.)**(5./6.)+ 

C  ((P4-l.)**2.+X**2.)**(5./6.)* 

C  ((P3-l.)**2.+Y**2.)**(5./6.)- 
C  2.*((P4+l.)**2.+X**2.)**(5./6.)* 
C  ((P3)**2.+Y**2.)**(5./6.)- 
C  2.*((P4-1.)**2.+X**2.)**(5./6.)* 
C  ((P3)*>t'2.+Y**2.)**(5./6.)- 
C  2.*((P4)**2.+X**2.)**(5./6.)* 

C  ((P3-l.)**2.+Y**2.)**(5./6.)- 
C  2.*((P4)**2.+X**2.)**(5./6.)* 

C  ((P3+l.)**2.+Y**2.)**(5./6.)+ 

C  4.*(P4**2.+X**2. )**(5./6.)* 

C  (P4**2.+Y**2,)**(5./6.)+ 

C  ((P4+1.)**2.+X**2.)**(5./6.)>K 
C  ((P4+l.)**2.+Y**2.)**(S./6.)+ 

C  C(P4-l.)**2.+X**2.)**(5./6.)* 

C  C(P4+l.)**2.+Y**2.)**(5./6.)+ 

C  ((P4+l.)**2.+X**2.)**<5./6.)* 

C  ((P4-l.)**2.+Y**2.)**(5./6.)+ 

C  ((P4-l.)**2.+X**2.)**(5./6.)* 

C  ((P4-l.)**2.+Y**2.)**(5./6.)- 
C  2.*((P4+l.)**2.+X**2.)**(5./6.)* 
C  ((P4)**2.+Y**2.)**(5./6.)- 
C  2.*((P4-l.)**2.+X**2.)**(5./6.)* 
C  ((P4)**2.+Y**2.)**(5./6.)- 
C  2.*((P4)**2.+X**2.)**(5./6.)* 

C  ((P4-l.)#*2.+Y**2.)**(5./6.)- 
C  2.*((P4)**2.+X**2.)**(5./6.)* 

C  ((P4+l.)**2.+Y**2.)**(5./6.) 

C  ) 

RETURN 

END 

REAL  FUNCTION  WTFUNCB(X) 

COMMON  /PARAMS/P3,P2,P3,P4 
REAL  P1,P2,P3,P4 
WTFUNCB=-(1.-ABS(X))*( 

C  2.*(Pl**2.+X**2.)**(S./6.)“ 
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C  ((Pl-l.)**2.+X**2.)**(5./6.)- 

C  ((Pl+l.)**2.+X**2.)**(5./6.) 

C  ) 

RETURN 

END 

:|c  >|c  :|c  ^  I):  *  sic  ><c  ^ i|c  >|t  4=  ****  if:  4^  **’)<’<<**  ^  * ’I' lie  *  4c  >1^  >1:  **:<<*  * 

REAL  FUNCTION  WTFUNCC(X) 

COMMON  /PARAMS/P1,P2,P3,P4 
REAL  P1,P2,P3,P4 
WTFUNCC=- ( 1 . -ABS (X) ) * ( 

C  2.*(P2*=|c2.+X**2.)**(5./6.)- 

C  ((P2-l.)**2.+X**2.)**(5./6.)- 

C  ((P2+l.)*>tc2.+X**2.)**(5./6,) 

C  ) 

RETURN 

END 

It:  :ic  4: 4c  lie  4c  *  *  4  i|c  4  )|cilc  lie  4cili  lie  lie  He  lit  lie  lie  lie  *  *  *  *  *  4s  «*  *  *  %  4=  4c  lie « i|t  i|c  :ic  4: 4c  ic  i|c  lie  *  %  *  « sic 

REAL  FUNCTION  WTFUNCD(X) 

COMMON  /PARAMS/P1,P2,P3,P4 
REAL  P1,P2,P3,P4 
WTFUNCD=-  ( 1 .  -  ABS  (X)  )  4C  ( 

C  2.4c((P3)4s4c2.+X4s4c2.)4c4s(5./6.)- 

C  ((P3-1.)4c4c2.+X4s4c2.)4c4.(5./6.)- 

C  ((P3+l.)ste4c2.+X4c4c2.)4c4<(5./6.) 

c  ) 

RETURN 

END 

4e4-4e4e4e4e4c4e4c4e4csieste4c4e4e4e4c4e4c4e4c4e4e4e4e4e4e4csie4e4e4c4e4s4s4c4e4e4e4c4e4c4e4estc4e4c4c4e4e 

REAL  FUNCTION  WTFUNCE(X) 

COMMON  /PARAHS/P1,P2,P3,P4 
REAL  P1,P2,P3,P4 
WTFUNCE=-  ( 1 .  -ABS(X) )  =k  ( 

C  2.sC((P4)4c4c2.+X4csK2.)4c4c(5./6.)- 

C  ((P4-l.)4c4e2.+X4cs|c2.)4:4e(5./6,)- 

C  (  (P4+1.)  4:4=2. +X4c4e2.)**  (5-/6.) 

C  ) 

RETURN 

END 

sfe  9|c)|(  j|c){e  sfcfc  9fe  9k  9f( 3fc  3f(  9fe  9|c  )|(  jfc  4e  4(:te ]fc  94c  V' 

SUBROUTINE  TRAPZD(FUNC,A,B,S,N) 

IF  (N.EQ.l)  THEN 

S=0.54<(B-A)4e(FUNC(A)+FUNC(B)) 

IT=1 

ELSE 

TNM=IT 
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DEL=(B-A)/TNM 
X=A+0.5*DEL 
SUM=0 . 

DO  200  J=1,IT 
SUM=SUM+FUNC(X) 

X=X+DEL 

200  CONTINUE 

S=0 . 5* (S+ (B-A) *SUM/TNM) 

IT=2*IT 

ENDIF 

RETURN 

END 

*  i|<  >11  !tc  4:  «  41  ***  >l<  Ik  !)<  **  <1 41  *  41  *  >|t  *  ifi  It:  *  >tc  *  lit ’ll  *  >)'****  >*< 

SUBROUTINE  qSIMP(FUNC,A,B,S) 

PARAMETER  (EPS=1.E-S,  JHAX=20) 

0ST=-1.E30 
0S=  -1.E30 
DO  300  J=1,JMAX 

CALL  TRAP2D(FUNC,A,B,ST,J) 

S=(4.4iST-0ST)/3. 

IF  (abs(s-os).lt.eps4:abs(os))  return 

OS=S 
OST=ST 
300  CONTINUE 

PAUSE  ’Too  many  steps.’ 

END 

4i4t4<4<4i4i4i4‘4i4<4<4i4i4i4i4i4<4i4<4i4‘4<4i4‘4i4i4i4i4i4t4i4i4i4i4<4i4i4:4i4i4<4t«4i4c4i4i’^4i4i4c 

SUBROUTINE  TRAPZDX(FUNC,A,B,S,N) 

IF  (N.EQ.l)  THEN 
S=0 . 5* (B-A) * (FUNC (A) +FUNC (B) ) 

IT=1 

ELSE 

TNH=IT 

DEL=(B-A)/TNH 

X=A+0.5*DEL 

SUH=0, 

DO  200  J=1,IT 
SUH=SUH+FUNC(X) 

X=X+DEL 
200  CONTINUE 

S=0 . 541  (S+  (B-A)  4CSUM/TNH) 

IT=2*IT 

ENDIF 

RETURN 

END 
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9(e  :|c  :fc  3|c  sfc  jfc  )|e  jtc  ^  :|t  %  )|c  3)c  sic  :4c  %  sfc  %  3k 

SUBROUTINE  QSIMPX(FUNC,A,B.S) 

PARAMETER  (EPS=l.E-5,  JMAX=20) 

0ST=-1.E30 
OS=  -1.E30 
DO  300  J=1,JMAX 

CALL  TRAPZDX(FUNC,A,B.ST,J) 
S=(4.*ST-0ST)/3. 

IF  (ABS(S-OS)  .LT.EPS-^ABS(OS))  RETURN 
OS=S 
OST=ST 
300  CONTINUE 

PAUSE  ’Too  many  steps.’ 

END 

SUBROUTINE  TRAPZDY(FUNC,A,B,S,N) 

IF  (N.EQ.l)  THEN 
S=0 . 5* (B-A) *  CFUNC(A)+FUNC(B) ) 

IT=1 

ELSE 

TNM=IT 

DEL=(B-A)/TNM 

X=A+0.5*DEL 

SUM=0. 

DO  200  J=1,IT 
SUM=SUM+FUNC(X) 

X=X+DEL 

200  CONTINUE 

S=0 . 5* (S+ (B-A) *SUM/TNH) 

IT=2*IT 

ENDIF 

RETURN 

END 

SUBROUTINE  QSIMPY(FUNC,A,B,S) 

PARAMETER  (EPS=l.E-5,  JMAX=20) 

0ST=-1.E30 
0S=  -1.E30 
DO  300  J=1,JMAX 

CALL  TRAPZDY(FUNC,A,B,ST,J) 
S=(4.*ST-0ST)/3. 

IF  CABS(S-OS).LT.EPS*ABS(OS))  RETURN 
OS=S 
OST=ST 
300  CONTINUE 
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PAUSE  ’Too  many  stops.’ 
END 
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